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Abstract 

This  paper  extends  past  analysis  of  an  optimal  source  distribution  around  a 
homogeneous  sphere  of  muscle  tissue  by  using  a  3-D  finite  difference  time  domain  (FDTD) 
scenario  in  which  an  anatomically  correct  human  head  model  is  irradiated.  It  first 
duplicates  the  analytical  solution  within  an  FDTD  space  using  an  FDTD  computer  code 
developed  at  Penn  State  University.  This  duplication  uses  a  9.45  cm  radius  sphere 
represented  in  an  FDTD  space  of  2.55  mm  cubic  cells.  FDTD  simulations  are  then 
peformed  on  four,  three,  and  two  layer  laminated  spheres,  designed  to  provide  simple 
models  of  a  head.  Finally,  four  simulations  were  performed  in  FDTD  on  the  human  head 
model  developed  at  Penn  State  from  an  MRI  scan  of  an  actual  human  head.  The 
comparison  of  analytic  simulations  to  the  FDTD  simulations  on  a  homogeneous  sphere 
showed  a  pixel  by  pixel  average  of  5.34%  error  between  the  two  with  a  standard  deviation 
of  7.84%.  The  layered  sphere  models  showed  considerable  spiking  at  the  two  poles  along 
with  a  small  amount  of  spiking  due  to  the  stair-step  approximation  of  the  spheres.  None  of 
these  spikes  increased  the  power  beyond  that  at  the  surface  and  hence  were  not  critical. 

The  simulations  on  a  true  human  head  showed  improvement  in  depth  due  to  the  low-loss  of 
the  bone  tissue.  This  study  demonstrates  that  microwave  hyperthermia  with  good 
resolution  is  possible  in  an  anatomically  correct  head  model. 
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1.  introduction 


Cancer  has  been  a  leading  cause  of  death  in  the  United  States  for  several  decades 
and  remains  so  today;  therefore,  it  is  a  leading  topic  of  research.  Current  researchers 
employ  several  methods  to  destroy  and  limit  the  growth  of  cancerous  tissue.  However,  all 
methods  contain  a  similar  characteristic;  they  destroy  the  healthy  tissue  as  well  as  the 
tumor.  One  promising  form  of  treatment  for  cancer,  for  limiting  the  damage  to  healthy 
tissue,  is  non-invasive  hyperthermia.  In  this  method  of  therapy  either  ultrasound  or 
electromagnetic  waves  are  inserted  into  the  tissue  from  many  different  locations.  These 
waves  induce  molecular  vibrations  that  heat  and  kill  the  cancer.  The  waves  pass  through 
the  tissue  and  cross  at  one  point  where  constructive  interference  occurs.  At  this  location, 
the  waveform  amplitude  is  significantly  higher  than  at  any  other  point  in  the  tissue.  Since 
the  amplitude  of  the  wave  dictates  the  amount  of  molecular  vibration  and  the  amount  of 
vibration  m  turn  dictates  the  amount  of  heat,  the  amplitude  controls  the  heat  produced.  In 
other  words,  the  location  of  constructive  interference,  with  its  increased  amplitude,  is 
heated  more  than  other  locations  in  the  tissue. 

The  current  drawback  to  non-invasive  hyperthermia  is  achieving  the  optimal 
solution  of  resolution  and  deep  penetration.  The  two  waveforms,  ultrasound  and 
electromagnetic  waves,  used  in  this  treatment  both  have  their  limitations.  Ultrasound 
provides  good  resolution,  however,  it  has  large  reflections  at  the  boundaries  of  the  different 
materials.  These  reflections  make  it  difficult  to  treat  tumors  near  dense  bone  tissue  or  open 
air  cavities,  such  as  tumors  in  the  head  or  chest.  Electromagnetic  waves,  on  the  other 
hand,  can  penetrate  these  substances  with  little  reflection,  but  high  frequencies  have  not 
been  able  to  penetrate  deeply  into  high  water-content  tissue  such  as  muscle. 

Problem  Statement 

This  research  focuses  on  using  high-frequency  electromagnetic  waves  as  part  of 
non-invasive  hyperthermia  in  treating  brain  cancers.  Although  deep  penetration  through  the 
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muscular  tissue  of  the  brain  has  been  impossible  in  the  past,  this  study  will  show  that  with 
the  combination  of  high  frequencies  and  constructive  interference,  resolution  at  depth 
improves  and  electromagnetic  waves  become  the  treatment  of  choice. 

Summary  of  Chapters 

This  thesis  is  broken  into  six  main  chapters,  introduction,  background,  theory, 
methodology,  results  and  conclusions.  The  introduction  you  have  just  completed  reading 
and  requires  little  explanation.  The  background  section  describes  the  work  that  has  been 
done  in  the  past  on  sphere  analysis,  FDTD  development,  and  material  characteristics  of 
biological  tissue.  The  theory  section  provides  a  brief  survey  of  the  theories  necessary  to 
complete  the  research  set  out  in  this  thesis.  This  includes  a  description  of  Rappaport’s 
optimization  process  and  Luebbers’  FDTD  code.  The  Methodology  section  describes  the 
modifications  made  on  Luebbers’  code  as  well  as  the  configuration  of  each  of  the  FDTD 
simulations.  The  results  section  describes  the  power  patterns  created  from  each  FDTD 
simulation.  Finally,  the  conclusion  section  provides  the  conclusions  which  can  be  drawn 
from  this  research  as  well  as  some  suggestions  for  future  work. 
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2.  Background 


A  significant  amount  of  research  has  evolved  the  tools  necessary  to  carryout  a  task 
as  difficult  as  analyzing  an  inhomogeneous  3-D  model  of  an  actual  human  head.  The  most 
important  is  the  work  by  Rappaport  to  derive  an  optimal  field  distribution  for  irradiating  a 
tumor  within  a  homogeneous  sphere  of  muscle  tissue.  This  research  demonstrates  the 
potential  of  microwave  hyperthermia.  To  apply  this  theory  to  the  inhomogeneous  structure 
of  a  human  head,  this  thesis  will  implement  the  finite  difference  time  domain  technique  for 
electromagnetics.  This  theory  approximates  Maxwell's  derivative  equations  as  finite 
differences  and  directly  solves  them  across  time  and  space.  Finally,  a  good  understanding 
of  the  electrical  properties  of  biological  tissues  was  necessary  to  model  adequately  the 
different  tissue  types  of  a  human  head. 

2.1  The  Sphere  Analysis 

Current  research  analyzes  only  perfect  spheres  of  homogeneous  tissue  [1, 2,3,4]. 

Dr.  Carey  Rappaport  at  Northeastern  University  leads  the  work  in  spherical  source 
radiation  for  hyperthermia  treatment  of  cancer.  In  his  work  he  has  dealt  solely  with  analytic 
solutions  to  homogeneous  spheres.  He  began  with  the  optimization  of  a  spherical  source 
distribution  to  irradiate  a  tumor  at  the  center  of  a  sphere  of  homogeneous  tissue.  Following 
the  analysis  at  the  center  of  a  sphere,  he  continued  on  to  investigate  heating  off-center 
targets  with  both  a  longitudinally  and  azimuthally  polarized  E-field. 

The  earliest  work  comes  from  Dr.  Rappaport  and  Dr.  Morgenthaler,  a  professor  at 
Massachusetts  Institute  of  Technology  [1].  They  place  a  tumor  at  the  center  of  a  sphere  of 
muscle  tissue  and  determine  the  maximum  radius  through  which  the  tumor  can  be  heated. 
By  determining  the  maximum  radius  with  an  ideal  distribution,  this  research  produces  a 
benchmark  other  methods  such  as  adaptive  nulling  can  try  to  obtain  but  never  improve.  To 
define  this  maximum  radius,  they  first  set  up  a  uniform  amplitude  and  longitudinally 
polarized  current  distribution  around  the  sphere.  They  then  calculate  the  radiation  pattern 
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within  the  sphere.  To  kill  only  the  cancer  cells  it  is  imperative  the  power  at  the  focus  be 
greater  than  at  all  other  points  within  the  sphere.  An  analysis  of  the  power  distribution 
which  the  surface  current  produces  reveals  that  the  point  where  power  at  the  surface  is  the 
same  as  that  at  the  focus  defines  the  maximum  depth  obtainable.  A  power  profile  along  the 
axis  of  the  sphere  is  computed  for  several  frequencies  from  100  MHz  to  2  GHz.  The  large 
range  of  frequencies  tested  was  to  examine  the  differing  characteristics  of  high  and  low 
frequencies.  Higher  frequencies  have  better  precision  but  exceptionally  small  penetration 
depths,  whereas  lower  frequencies  have  excellent  penetration  but  poor  precision.  This 
analysis  shows  the  optimal  frequency  to  be  915  MHz.  This  occurs  because,  due  to  the 
electrical  characteristics  of  the  muscle  tissue,  915  MHz  actually  has  better  precision  and 
better  depth  than  the  next  frequency  down,  700  MHz.  Rappaport  and  Morgenthaler 
compute  the  maximum  radius  to  be  just  under  9  cm  for  915  MHz. 

To  increase  this  maximum  radius  obtainable,  Rappaport  and  Morgenthaler  turn  to  a 
modal  analysis  in  which  they  add  higher  order  spherical  harmonics  to  the  field  distribution 
around  the  sphere.  The  modal  analysis  allows  higher  order  modes,  which  have  no  power 
at  the  center,  to  destructively  interfere  near  the  edges  around  the  equator  of  the  sphere. 
Since  surface  power  is  highest  at  the  equator  of  the  sphere,  this  technique  attempts  to  force 
the  surface  power  around  the  equator  to  decrease,  spreading  that  energy  toward  the  poles. 
However,  even  with  these  higher  harmonics,  the  theoretical  maximum  radius  through 
which  a  tumor  can  be  heated  increased  to  only  9.45  cm  from  9  cm  at  915  MHz.  Although 
this  is  close  to  a  typical  head  size,  it  assumes  an  ideal  source;  thus,  achieving  it  in  practice 
would  be  impossible. 

Rappaport,  with  his  graduate  student  Pereira,  then  turned  toward  radiating  off 
center  targets.  The  one  difference  between  the  off-center  heating  scheme  and  the  center 
heating  configuration  lies  in  polarization.  For  a  tumor  located  at  the  center  of  a  sphere,  the 
polarization  of  the  field  at  the  focus  is  not  a  factor.  The  polarization  of  the  field  distribution 
on  the  surface  must  be  in  the  same  direction  over  the  entire  surface;  however,  the  direction 
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of  this  polarization  is  unimportant.  If  there  is  the  desire  to  have  another  polarization  at  the 
center  of  the  sphere,  the  source  distribution  need  only  to  be  rotated  about  the  outside  of  the 
sphere.  If  the  field  distribution  rotates  in  the  case  of  an  off-center  focus,  the  location  of  the 
focus  would  move. 

The  off-center  focus  requires  two  separate  derivations,  one  for  each  polarization. 
To  define  these  polarizations,  it  is  necessary  to  define  two  axes  of  the  sphere.  The  first 
axis  defined  passes  directiy  through  the  center  of  the  tumor  and  is  therefore  called  the 
centered  axis.  The  second  axis  simulations  perpendicular  to  the  first  and  is  labeled  the  off- 
center  axis.  Figure  1.1  illustrates  the  two  axes’  definitions.  The  polarization  is  then  set 
parallel  to  the  off-center  axis  (azimuthal  polarization)  or  parallel  to  the  centered  axis 
(longitudinal  polarization). 


Figure  2.1:  This  figure  illustrates  the  definition  of  the  centered  and 

off-center  axes. 

Rappaport  and  Pereira's  paper  looks  at  longitudinal  polarization.  They  closely 
followed  what  Rappaport  had  done  with  Morgenthaler  [2,3].  The  difference  is  the  location 
of  the  origin.  For  the  off-center  study,  they  locate  the  origin  of  analysis  at  the  tumor 
location.  Computing  the  fields  that  correspond  to  a  source  at  that  location,  they  set  up  a 
field  distribution  that  focused  at  the  off-center  location.  In  other  words,  they  place  the 
center  of  the  spherical  wave  distribution  at  the  tumor  location  and  determine  the 
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corresponding  field  at  locations  on  the  surface  of  the  sphere.  With  a  modal  analysis  similar 
to  the  one  used  for  heating  a  centered  target,  they  are  able  to  increase  precision  at  the  off 
center  location.  For  heating  a  location  at  six-tenths  of  the  radius,  they  are  able  to  penetrate 
through  a  12  cm  radius  sphere  at  915  MHz. 

Upon  completing  the  longitudinal  analysis  with  Pereira,  Rappaport  published  a 
paper  looking  at  the  azimuthal  polarization  of  the  off-center  target[4].  The  azimuthal 
polarization  analysis  entails  a  similar  analysis  as  the  longitudinal,  except  polarization  of  the 
field  is  parallel  to  the  off-center  axis.  This  polarization  is  not  as  effective  as  longimdinal 
polarization.  For  treatment  at  six-tenths  of  the  radius  Rappaport  is  only  able  to  penetrate 
through  a  9  cm  radius  sphere  at  915  MHz.  It  is  possible  to  understand  the  difference  from 
an  intuitive  sense.  The  maximum  power  addition  comes  from  the  sources  near  the  equator 
of  the  source  distribution,  for  these  are  the  strongest  sources.  For  the  longitudinal  case, 
these  sources  all  lay  equal  distances  from  the  tumor.  In  azimuthal  polarization  some  of 
these  tangential  sources  are  located  further  than  the  sphere's  radius  away.  Transmission 
through  lossy  tissue  attenuates  the  field  significantly,  reducing  the  power  contribution  of 
these  sources,  which  in  turn  reduces  the  maximum  radius  through  which  tumor  treatment  is 
possible. 

2.2  FDTD  Background 

The  computational  finite  difference  time  domain  (FDTD)  method  also  provides  a 
method  for  analysis  of  electromagnetic  wave  propagation  through  tissue.  FDTD  breaks 
time  and  space  into  finite  increments.  It  then  applies  Maxwell's  equations  directly  to  this 
grid,  with  the  derivatives  approximated  as  differences.  Dr.  Raymond  Luebbers  of  Penn 
State  University  developed  a  general  purpose  3-D  FDTD  code  that  is  extremely  versatile 
and  easy  to  use  [5].  It  includes  a  graphical  user  interface  (GUI)  for  inputting  the 
geometries  as  well  as  displaying  results.  Although  developed  for  backscatter  predictions 
and  antenna  design  problems,  this  codes  versatility  allows  straightforward  modification  to 
handle  a  biological  tumor  treatment  scenario.  In  addition,  Mr.  David  Steich  has  produced  a 
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human  head  model  capable  of  use  with  this  FDTD  code.  This  geometry  mesh  was  created 
from  an  MRI  scan  of  an  actual  human  head.  It  is  important  to  note  that  this  is  an 
approximate  model  but  will  prove  to  be  adequate  forjudging  the  effects  of  the 
inhomogeneities  of  the  human  head. 

Om  Ghandi  and  Carl  Dumey  at  the  University  of  Utah  have  done  a  great  deal  of 
research  in  the  field  of  biological  FDTD  research.  Their  research  concerns  the  biological 
effects  of  electromagnetic  waves,  rather  than  the  treatment  of  biological  tumors.  However, 
the  methods  they  use  to  predict  field  propagation  through  tissue  directly  apply  to  the 
treatment  scenario  [6]. 

The  research  done  by  Ghandi  and  Dumey  has  two  important  lessons  that  pertain  to 
my  research.  First,  It  demonstrates  the  versatility  of  the  FDTD  code  as  a  tool  for  analyzing 
electromagnetic  effects  in  biological  tissue.  The  grid  set  up  in  FDTD  may  contain  several 
different  electric  or  magnetic  materials  in  virtually  any  configuration.  This  allows  relatively 
easy  modeling  of  the  complex,  inhomogeneous  structure  of  the  human  body.  The  study 
also  demonstrates  the  limits  of  actual  applicators.  Rappaport's  studies,  as  well  as  mine, 
assume  the  ability  to  produce  an  E-field  distribution  exactly  as  desired.  The  study  by 
Dumey  et.  al.  shows  that  in  practice  this  wiU  be  nearly  impossible  because  of  the  capacitive 
effects  that  occur  due  to  the  high  conductivity  of  the  muscle  tissue.  The  muscle  acts  as  a 
second  plate  of  a  parallel  plate  capacitor,  hence  the  field  does  not  behave  as  it  would  for  a 
simple  radiating  dipole  [6].  It  is  essential  to  keep  this  drawback  in  mind  when  doing 
research  with  ideal  distributions.  The  ideal  distribution  will  never  be  obtainable;  therefore, 
to  have  hope  of  a  workable  system,  the  ideal  maximum  depth  must  be  somewhat  deeper 
than  real  applications  will  ever  encounter.  This  will  allow  for  the  loss  of  depth  that  will 
come  as  a  result  of  the  inability  to  create  a  source  distribution  exactly  as  predicted  in  the 
ideal  scenario. 
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2.3  Electrical  Properties  of  Human  Tissue 

Curtis  Johnson  and  Arthur  Guy  published  a  paper  in  1972  from  which  most 
research  in  1993  has  obtained  model  electrical  characteristics  for  human  tissue  [7].  In  this 
paper,  they  "discuss  the  electric  properties  of  differing  forms  of  tissue.  The  study  shows 
that  muscle  tissue  is  the  lossiest  tissue  of  the  human  body.  That  is,  electromagnetic  waves 
attenuate  more  quickly  through  muscle  than  any  other  form  of  tissue.  This  characteristic 
dictates  muscle  tissue  as  the  tissue  of  choice  for  calculating  the  performance  of 
electromagnetic  treatment.  This  provides  a  worst-  case  scenario.  In  a  real  head  model, 
there  is  relatively  low-loss  bone  tissue  and  open  air  cavities.  This  means  more  boundary 
reflections,  but  also  means  less  loss  as  the  wave  propagates  through  the  tissue. 

In  1980  Stuchly  and  Stuchly  published  a  more  complete  set  of  dielectric  properties 
of  biological  tissues  [8].  This  list  breaks  down  the  material  characteristics  more 
completely,  giving  the  permitivities  of  bone,  muscle,  white  matter,  gray  matter,  etc.  This 
indicates  the  brain  itself  also  has  lower  loss  than  muscle  tissue.  Thus,  there  remains  the 
possibility  that  a  small  deep-set  tumor  within  an  actual  human  head  model  may  be  treatable, 
given  an  optimal  source  distribution. 

The  research  done  to  date  indicates  that  although  there  are  certain  limitations  to 
electromagnetic  hyperthermia  treatment  for  cancer  reduction,  the  technique  still  holds  many 
benefits  over  other  forms  of  hyperthermia  and  is  worthy  of  continued  study.  The  limited 
precision  at  deep  locations  is  a  shortcoming  of  electromagnetic  waves  that  is  avoidable,  as 
seen  in  Rappaport's  work  and  will  be  shown  here.  In  any  event,  it  is  an  area  that  clearly 
deserves  more  research.  The  inability  of  ultrasound  to  penetrate  bone  and  air  pockets 
makes  it  impossible  to  use  ultrasound  waves  for  treatment  in  the  skull  and  chest  region. 

The  low  loss  of  bone  tissue  and  air  cavities  leaves  electromagnetics  a  strong  contender  for 
future  research  in  cranial  and  cardial  tumors.  Optimization  of  the  source  distribution  using 
an  FDTD  analysis  of  an  inhomogeneous  head  can  be  done  to  develop  a  microwave  system 
capable  of  heating  only  the  cancerous  tissue,  leaving  the  remaining  tissues  unharmed. 
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3.  Theory: 


There  are  two  main  theories  used  throughout  this  research.  The  first  is  the  Finite 
Difference  Time  Domain  (FDTD)  technique  for  electromagnetics.  This  theory  uses  finite 
difference  approximations  for  the  derivatives  of  Maxwell's  equations.  It  then  cycles 
through  time  steps  and  solves  these  equations  at  every  location  within  the  FDTD  space.  In 
this  way  predicting  the  propagation  of  electromagnetic  waves  through  complex  media  is  a 
matter  of  solving  millions  of  simple  calculations.  The  second  theory  to  be  understood  is 
the  solution  to  the  vector  wave  equation  in  spherical  coordinates  and  how  it  applies  to  the 
optimization  technique  developed  by  Rappaport.  This  theory  defines  the  optimal  source 
distribution  for  irradiating  deep  set  tumors. 

3.1  Finite  Difference  Time  Domain  Technique 

Unlike  Rappaport's  work,  this  thesis  will  not  use  an  analytic  approach  for  solving 
for  the  electromagnetic  fields.  Instead,  it  uses  a  numerical  approximation  technique  known 
as  the  Finite  Difference  Time  Domain  (FDTD).  FDTD  approximates  the  differential  form 
of  Maxwell's  equations  with  finite  differences.  To  do  this,  it  is  necessary  to  break  both 
time  and  space  down  into  small  increments.  The  method  then  directly  applies  Maxwell's 
equations  to  each  of  these  cubes  at  each  time  increment  thereby  propagating  the 
electromagnetic  fields  through  the  space.  This  allows  application  of  any  material 
characteristic  to  any  cube  within  the  space,  making  analysis  of  inhomogeneous  bodies  in 
the  near  field  exceptionally  simple.  The  following  derivation  is  from  [9]. 

As  stated  above,  FDTD  is  simply  a  finite  difference  application  of  Maxwell's 
equations  in  differential  form.  Let  us  therefore  begin  at  Maxwell's  equations  in  differential 
form  in  the  time  domain. 
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Where 


V  X  E  =  /  dt 

WxR  =  dD/ dt  +  J 
VD  =  p 
VB  =  0 


D  =  e*E 
B  =  p*H 

It  is  important  to  stop  and  point  out  that  the  constitutive  relations  here  are 
convolutions  and  not  simple  multiplications.  This  is  because  these  equations  apply  to  the 
time  domain,  and  not  the  frequency  domain.  Because  the  equations  in  the  frequency 
domain  are  multiplications,  in  the  time  domain  they  become  convolutions.  If  we  deal  with 
only  one  frequency  and  assume  a  scalar  permitivity  and  permeabihty,  multiphcations 
replace  these  convolutions. 

Before  going  any  further,  we  can  reduce  the  derivation  to  only  the  two  curl 
equations  by  noting  that  the  divergence  equations  are  redundant.  To  show  these  are 
unnecessary,  start  by  taking  the  divergence  of  each  of  the  curl  equations. 

V-(VxE)  =  V-  ~  (1) 

V  ut  J 

V-(VxH)  =  V-[^^  +  jj  (2) 


Now  by  using  V  •  V  x  A  =  0  it  is  easy  to  see  eq(l)  and  eq(2)  transform  into 

,3) 

0^3(V  D)^y  ^  (4) 

dt 


■  +  V-J  (4) 


Eq  (3)  implies 


V  •  B  =  constant  wrt  time  (5) 


By  using  the  continuity  equation,  V  •  J  +  ^  =  0,  eq(4)  can  be  rewritten  as 

ot 

—[(V  •  D)  -  p]  =  0  which  implies 

at 

V  •  D  -  p  =  constant  wrt  time .  (6) 
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Since  there  will  be  no  sources  or  currents  at  t=0,  V  •  B  =  0  and  V  •  D  -  p  =  0  at  initial  time. 
Eq(5)  and  eq(6)  show  that  these  equalities  must  hold  for  all  time.  Since  we  have  just 
derived  the  divergence  equations  from  the  curl  equations,  the  divergence  equations  are 
redundant. 

Now,  to  allow  for  magnetic  loss  we  will  introduce  a  magnetic  current  term 
M  =  ct‘H,  then  rearrange  the  first  two  equations  to  produce  two  differential  equations  of 
the  form 

^  =  --(VxE)-— H  (7) 
dt  11  fl 

■^  =  -(VxH)--E  (8) 

In  taking  this  step,  the  p  and  £  were  removed  from  within  the  time  derivative.  This 

requires  that  the  permittivity  and  permeability  cannot  be  functions  of  time  (or  frequency). 

Since  in  human  tissue  this  is  not  the  case,  the  analysis  will  use  only  one  frequency  with 

p  and  £  given  by  a  constant.  This  frequency  limitation  converts  the  convolution  in  the 

constitutive  relations  to  a  multiplication  with  a  scalar,  and  allows  the  removal  of  the 

permitivity  and  permeability  from  within  the  derivative. 

The  total  E-field  can  be  written  as  the  sum  of  an  incident  field  and  a  scattered  field. 

The  incident  field  is  defined  as  propagating  through  free  space  everywhere  throughout  the 

problem  space  and  the  scattered  field  is  defined  as  propagating  through  free  space  outside 

the  scatterer  and  propagating  through  the  media  inside  the  scatterer.  This  sum  is  then  be 

inserted  into  eq(7)  and  eq(8)  to  produce 

^^*L±Mi  =  _l(Vx(E''+EO)-— (H'+HO  (9) 
dt  p ''  ^  p 

and 

^(E  +E  )=l/vx(H‘'+ff))--(E'+EO  (10). 

dt  £ '  ’  £ 

Since  the  incident  field  is  defined  in  free  space  it  has  to  obey  the  equations 

and 
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(12). 


VxH' 

The  incident  fields  are  generally  known  quantities.  Therefore,  to  produce  an  equation  for 

the  scattered  field,  subtract  the  value  for  the  incident  field  from  eq(7)  and  eq(8).  At  the 

same  time,  the  corresponding  values  in  eq(ll)  and  eq(12)  replace  the  curl  of  the  incident 

field.  These  substitutions  produce  two  equations  for  the  scattered  fields  of  the  form 

—  =  +  -(VxHO  (13) 

dt  £  £  £  dt  £ 

and 

—  =  --(VxEO  (14). 

dt  II  II  n  dt  II 

These  differential  equations  are  then  put  into  finite  difference  form  by  implementing 
the  definition  of  a  derivative, 

f(x,t^)-fix.t^)  f{x,t^)- fjx.t^) 

dt  At  — >  0  At  At 

f{X^,t)-f{Xy,t)  fiX2,t)-f{X^,t) 

dx  Ax  -^0  Ax  Ax 

Note  that  the  del  operator  in  eq(13)  and  eq(14)  contains  the  spatial  derivative.  After 
producing  these  finite  difference  equations,  it  is  a  simple  matter  of  breaking  down  the 
equations  into  x,  y,  and  z  components  and  implementing  the  finite  difference  equations  on  a 
computer. 

Implementing  the  differential  form  of  the  time  derivative  on  eq(9)  and  eq(lO) 

produces  the  following  difference/derivative  equations. 

e(E^’'’  -E*’”-^)+<TAfE^'''  =-(TAtE'-"-(e-eJAtE‘'''  +At(VxH '""2)  (17) 

I 

’”■')+ CT*  Atff =-0r*AtH''''-(/i-;UjAfH''-At(VxE'^2)  (18) 

The  introduction  of  the  superscript  n  represents  the  nth  time.  In  addition  the  shorthand 
notation  E  represents  the  derivative  with  respect  to  time.  The  "leapfrog"  method  of 
updating  the  fields  is  also  evident  in  these  equations  in  the  half  time  step  offset  of  the 
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electric  and  magnetic  fields.  The  "leapfrog"  technique  will  be  discussed  more  fully  in  the 
discussion  of  implementing  the  code. 


Solving  eq(13)  and  eq(14)  for  the  scattered  fields  yields 

f  crAf 


E"’"  = 


£+(7Atj 

i£-£o)At 

£  +  (jAt 


iE<.«  + 


£+(jAt, 

'  At 
,£  +  crAr 


(VxH 


f  u  ^ 

_  h*-  Jj-f.n- 

+  G*  At ) 

(At 

H  +  G*At 

A  computer  code  then  solves  these  equations  across  time  and  space  to  predict 
electromagnetic  interactions  with  scattering  objects.  The  methodology  section  of  this  thesis 
discusses  implementation  of  these  equations  in  complete  detail. 

Before  implementation  of  FDTD,  the  cell  size  and  time  step  size  have  to  be 
determined.  In  general,  the  cell  size  determines  the  accuracy.  The  cells  have  to  be  at  least 
one  tenth  of  the  shortest  wavelength  within  the  materials.  The  smaller  the  cell  size,  the 
more  accurate  the  results  wiU  be.  If  the  scattering  object  has  ramped  or  spherical  surfaces, 
a  quantization  error  will  result.  This  error  is  commonly  referred  to  as  the  "staircase"  error. 
Smaller  cells  can  reduce  this  error  by  producing  a  more  accurate  rendering  of  the  actual 
geometry.  Reducing  the  cell  size  has  its  drawbacks,  however.  Reduction  of  the  cell  size 
requires  subsequent  reduction  of  the  time  step  size.  With  smaller  time  step  increments,  the 
code  must  cycle  through  more  time  steps  to  allow  for  the  same  amount  of  time  and  the  same 
field  propagation  distance. 

The  cell  size  determines  the  size  of  the  time  step  to  ensure  stability  in  the  solution. 
The  FDTD  simulation  can  only  propagate  the  field  across  one  cell  in  each  time  step.  Too 
long  of  a  time  step  forces  the  FDTD  simulation  to  propagate  the  field  across  multiple  cells 
in  one  time  step,  which  it  is  unable  to  do.  Thus,  the  solution  becomes  unstable.  In  order 
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to  determine  what  the  optimal  time  step  is,  we  turn  to  the  Courant  stability  condition  [9] 
which  states  that  for  a  propagation  speed  v  within  the  material,  the  time  step  must  be 


chosen  to  satisfy  the  condition. 


A  1  1  1 

[(Axf  (Ayf  {Azf 


This  condition  specifies  the  longest  time  step  that  will  maintain  stability.  The  use  of  shorter 
time  steps  is  possible,  but  show  little  improvement  and  are  not  worth  the  added 
computation  time  [9]. 

From  this  equation  it  can  be  seen  that  to  ensure  stability,  the  time  step  size  must  be 
reduced  when  the  ceU  size  is  reduced.  In  addition,  raising  the  frequency  forces  reduction 
of  the  cell  size  due  to  the  shortening  of  the  wavelengths;  thus,  reduction  of  the  time  step  is 
necessary.  To  produce  the  best  results  in  the  shortest  time,  it  is  necessary  to  weigh  the 
accuracy  desired  at  a  given  frequency  against  the  computation  time  necessary  to  set  up  the 
optimal  simulation  configuration.  The  methodology  section  describes  the  simulation 
configuration  which  this  research  uses. 


3.2  Source  Optimization 

The  source  optimization  developed  by  Rappaport  and  Morgenthaler  uses  the 
solution  to  the  vector  wave  equation  in  spherical  coordinates.  Appendix  A  gives  a  complete 


A  =  j„mp: (.cosd)  .  m(t>.  (23)  [10] 
mn  sin 


Where 


Before  going  further,  it  is  important  to  note  a  few  of  the  dependencies  within  this 
equation.  The  spherical  Bessel  function  of  order  n  predominately  controls  the  radial 


dependence.  The  m  and  n  modes  of  the  Legendre  function  dictate  the  theta  dependence. 
Finally,  the'  (j)  dependence  varies  as  sinm^.  Figure  3.1  describes  the  location  of  all  of 
these  components 


Figure  3.1:  The  Spherical  coordinate  system  as  used  in  this  thesis 

Rappaport  and  Morgenthaler  began  at  eq  (22)  to  find  the  optimal  source 
distribution.  The  following  description  of  this  optimization  comes  from  [1].  Solving 
eq(22)  for  the  derivatives  and  using  only  the  m=0  modes  produces 


16 


E(R)  ^-^fA^^^P^j.mp.icose) 


n=l 


+ft4 


kR 


h-m)-—jnikR) 

kR 


(24) 


*[n  cos  (cos  6)  -  nP„_^  (cos  d)]\  ^— 

vsin0 


Figure  3.2  shows  a  plot  of  the  power  across  a  vertical  cut  through  the  sphere.  This 
is  a  plot  of  the  power  with  only  the  first  mode  excited.  It  is  clear  from  this  plot  that  the 
power  is  concentrated  near  the  equator.  Increasing  the  power  near  the  poles  will  allow 
reduction  of  the  power  near  the  equator,  without  reducing  the  power  at  the  focus.  In  other 
words,  the  power  at  the  equator  needs  to  be  spread  out  toward  the  poles. 

Using  only  the  m=0  modes  spreads  the  field  as  evenly  around  the  sphere  as 
possible.  Higher  order  m  modes  would  force  a  sinusoidal  variation  in  the  Phi  direction. 
This  variation  would  produce  nulls  at  evenly  placed  intervals,  depending  on  the  value  of  m. 
These  nulls  would,  in  turn,  force  higher  power  at  the  peak  areas  on  the  surface  to  get  the 
same  amount  of  power  to  the  center  due  to  the  absence  of  power  at  the  nulls.  Thus,  unless 
there  was  for  some  reason  a  location  where  the  E-field  had  to  be  kept  at  zero,  it  would  be 
counterproductive  to  use  higher  order  m  modes. 

The  higher  n  modes,  however,  can  be  used  to  our  advantage.  Because  (kR) 
varies  as  (kR)"  for  small  values  of  R,  the  n=l  mode  is  the  only  mode  that  has  power  at 
r=0,  all  other  modes  are  zero  at  the  center.  Figure  3.3  shows  this  characteristic  of  higher 
order  n  modes.  It  is  a  plot  of  the  first  three  n  modes  peak  E-field  value  versus  radius, 
where  R=0  is  the  center  of  the  sphere  and  increasing  R  goes  out  along  an  equatorial  cut  of 
the  sphere.  Therefore,  addition  of  higher  order  n  modes  out  of  phase  with  the  first  mode, 
will  allow  their  use  to  reduce  the  power  on  the  surface  without  affecting  the  power  at  the 
center. 
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Power  Distribution  Across  Central  YZ  Plane  with  only  the  N=1  Mode. 


Power  (W/m^2) 


Z  Location 

Figure  3.2:  Plot  of  the  power  distribution  across  the  central  YZ  cut  with  only  the  N=1 

mode  excited. 
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E-Field  (V/m) 


Theta  Component  of  First  Three  Modes  Versus  Radius. 


Figure  3.3:  Peak  E-field  strength  versus  radius  for  the  first  three  n  modes.  R=0 

represents  the  center  of  the  sphere. 

However,  optimization  does  not  require  the  use  of  the  even  n  modes.  This  can  be 

seen  from  looking  at  the  instantaneous  E-field  values  along  the  surface  of  the  sphere.  For 
1 

large  radii,  the  —  dependence  reduces  the  r  component  and  the  6  component  becomes 

—  ^ 
the  dominating  factor.  Therefore,  to  reduce  the  power  levels  at  the  surface  (large  R)  the  d 

component  is  the  component  most  important  to  reduce.  Since  we  are  dealing  only  with  the 

m=0  mode,  there  is  no  6  dependence  in  the  (j)  direction,  see  eq(22).  Therefore,  the 

A 

optimization  must  deal  with  the  6  dependence  in  the  6  direction. 
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Instantaneous  E-Field  Stength  (V/m) 


N=2  Mode  Versus  Angle  Theta. 


Figure  3.4:  This  flgure  is  a  plot  of  the  instantaneous  E-field  versus  angle  theta 
for  the  n=2  mode  at  R=9.45  cm  (along  the  surface  of  the  sphere). 

Eq(19)  shows  the  6  dependence  to  be 

dPJsSSS.,  (25) 
dd 

P„  (cos  d)  is  an  even  function  of  6  for  even  values  of  n;  therefore,  its  derivative  will  be  an 

7t 

odd  function  of  6  and  will  have  a  null  at  ^  -  y-  Figure  3.4  demonstrates  this  in  a  plot  of 

A 

the  instantaneous  6  component  of  the  E-field  for  the  n=2  mode  at  a  radius  of  9.45  cm  (the 
surface  of  the  sphere).  Since  these  even  order  modes  have  a  null  at  the  equator,  where  we 
hope  to  reduce  the  power  level,  they  will  be  of  no  help  in  optimization.  In  addition,  since 
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Instantaneous  E-Field  Strength  (V/m) 


1.5 


N=3  Mode  Versus  Angle  Theta. 


Figure  3.5:  This  is  a  plot  of  the  instantaneous  E-field  versus  angle  theta  for  n=3 

mode  at  R=9.45  cm. 

they  are  odd  functions,  when  the  fields  subtract  in  the  lower  hemisphere  they  will  add  in 
the  upper  and  vice-versa. 

The  Legendre  functions  of  odd  order  n,  however  are  odd  functions  of  d. 

7t 

Therefore,  their  derivatives  are  even  functions  of  theta  and  will  have  a  peak  at  ^  =  "j- 

Figure  3.5  illustrates  this  idea  in  a  plot  of  the  instantaneous  E-field  for  the  n=3  mode  versus 
angle  theta  along  the  surface  of  the  sphere.  From  this  plot  it  can  be  seen  how  the  odd  order 
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n  modes  can  be  used  to  reduce  the  power  near  the  equator.  It  is  the  same  idea  as  Fourier 
analysis  of  periodic  time  signals  to  produce  a  pulse,  only  it  is  using  spatial  harmonics  rather 
than  time  harmonics.  The  reasons  given  above  demonstrate  the  optimization  should  deal 
solely  with'the  odd  order  harmonics  of  index  n.  It  is  interesting  that  because  the  r 
dependence  varies  with  P^{co&d),  its  even  order  harmonics  would  be  preferential  for 

optimization,  but  because  they  are  less  dominate,  we  deal  only  with  the  6  component. 

The  optimization  plan  was  to  add  these  first  three  odd  order  harmonics  in  such  a 
way  to  reduce  the  power  at  the  equator  as  much  as  possible.  To  do  this,  start  at  the  6 
dependence  of  the  derivative  of  the  Legendre  function  for  n=l,3  and  5. 

~  -sin0 

Eg. — — cos^  0sin0  +— sin0  (26) 

03  2  2 

- cos  0COS0H - cos  0sm» - smw 

8  4  8 

Using  trigonometric  identities,  rewrite  these  equations  in  terms  of  odd  multiples  of  6. 

Eg^  — sin0 

Egj  ~-ysin30  +  |sin0  (27) 

„  315  .  105  ...  30  .  . 

Ea< - Sin  50 - sin  30 - sin0 

128  128  128 

It  is  now  possible  to  rewrite  the  three  equations  (27)  into  one  equation  describing 
the  0  dependence  as 

/(0)  =  sin  0  +  5j  sin  30  +  02  sin  50  (28) 

Next,  solve  this  equation  for  the  B  coefficients  to  minimize  the  maximum  value  of  {f(d)f 
from  0  to  7t.  This  produces  coefficients  of  value  =.2365  and  fij  =.0640.  The  resulting 
form  of  the  equation  has  three  peaks  at  the  maximum  value  .685  of  the  previous  maximum. 
Figure  3.6  shows  a  plot  of  the  equation  for  the  first  three  odd  harmonics  as  well  as  for  the 
first  harmonic  alone. 
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Power  (W/m'^2) 


-  Single  Mode  Theta  (Degrees) 

—  First  Three  Modes 

Figure  3.6:  This  figure  is  a  plot  of  (/(0))^  for  the  first  three  odd  modes 
compared  to  the  plot  of  only  the  first  mode. 

From  here  it  can  be  seen  how  this  could  be  extended  out  in  terms  of  SIN’s  of 
higher  multiples  of  d  to  reduce  the  ripple  further.  However,  because  there  is  a  great  deal 
of  error  inherent  in  determining  the  material  characteristics  of  the  head,  this  decreased  ripple 
will  not  not  be  worth  the  extra  time  required  for  analysis. 
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The  next  step  is  to  add  in  the  spherical  Bessel  function  scahng  factors  on  each 

Tl 

mode.  This  is  done  by  introducing  the  coefficients  J„=\  j„_i  (kR) - (kR)  . 

L  kR  _ 

Combining  this  definition  with  the  B  coefficients,  eq  (28)  for  f{d) ,  and  the  three  equations 
(27)  produces  a  system  of  equations  for  the  /  coefficients  of  the  form 


U  --  +/<  - 


=  0.2365.  (29) 


=  0.0640 


The  solution  for  this  system  is  =  -0.9509,  =  -0.1 148,  and  /j  =  -0.0260. 


The  first  J  coefficient  is  reduced  from  unity  due  to  the  sin  9  terms  in  the  higher 
order  modes.  Thus,  the  power  at  the  center  of  the  sphere  is  reduced  by  the  same  amount 
(the  higher  order  modes  due  not  contribute  to  the  power  at  the  center).  Division  of  these 

f 

three  coefficients  by  -0.9509  produces  a  coefficient  equal  to  unity  and 

^3  =  0.1207,  and  =  0.0288.  Thus,  the  power  at  the  surface  will  be  72%  the  value  of  a 


single  mode  excitation. 

From  here  the  optimization  changes  based  on  the  characteristics  of  the  material  and 


the  frequency  used.  This  research  uses  a  frequency  of  915  MHz  which  corresponds  to  a 
relative  permittivity  in  muscle  tissue  equal  to  5 1  and  a  conductivity  of  1.26.  These  values 
are  used  to  solve  for  the  value  of  k  and  substituted  back  into  the  definition  of  the  J 


coefficients  to  find  values  for  the  A  coefficients.  In  solving  for  the  optimal  value  of  the  A 
coefficients,  the  equations  change  for  each  choice  of  radius  R,  thus  the  optimization  is  an 
iterative  process.  When  these  iterations  are  complete,  the  optimal  A  coefficients  for  a 
homogeneous  sphere  of  muscle  tissue  are  found  to  be  Al=1.21  +  0.89j,  A3=-.19  -  0.05i, 
and  A5=.05  -  .02i. 


24 


4.  Methodology 


This  research  entails  setting  up  three  distinct  simulations.  The  first  is  an  analytic 
simulation  that  serves  as  a  comparison  for  the  FDTD  results  to  validate  the  FDTD  code. 

For  this  simulation,  a  FORTRAN  program  solves  the  spherical  vector  wave  equation  on  a 
homogeneous  sphere  of  muscle  tissue.  A  series  of  FDTD  simulations  on  laminated  spheres 
comprise  the  second  set  of  simulation.  This  set  provides  insight  into  the  effects  of 
inhomogeneities  on  the  propagation  of  the  ideal  source  radiation.  A  series  of  FDTD 
simulations  on  an  actual  model  of  a  human  head  developed  from  an  MRI  scan  at  Penn  State 
comprise  the  final  set  of  simulations. 

4.1  Spherical  Vector  Wave  Equation  Computer  Code 

The  first  step  entails  duplicating  Rappaport's  analytic  results  to  use  as  a  control  with 
which  to  compare  FDTD  solutions.  A  FORTRAN  code  solves  Rappaport's  optimized 
solution  in  a  space  that  has  the  same  grid  size  as  the  proposed  FDTD  space.  Solving  this 
equation  across  a  grid  identical  to  the  one  used  in  the  FDTD  simulation  allows  easy 
comparison.  Figure  4.1  shows  the  flow  diagram  of  the  code. 
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Figure  4.1:  This  flgure  shows  the  flow  diagram  of  the  FORTRAN 
code  to  compute  the  analytic  solution  to  spherical  vector  wave 

equation. 
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The  main  body  of  the  code  steps  through  x,  y,  and  z  increments  each  running  from 
one  to  seventy-five.  As  in  the  FDTD  grid  that  will  be  used,  each  'cell'  is  cubic  with  a 
length  of  .002554  m.  As  the  code  steps  through  each  cube  in  the  grid,  it  computes  the 
angle  6  down  from  the  z  axis  and  the  angle  (j)  off  the  x  axis,  as  well  as  the  distance,  R, 
from  the  center  of  the  cube  of  interest  to  the  center  of  the  sphere,  that  is  located  at  x=37.5, 
y=37.5,  and  z=38.  The  code  then  calls  the  subroutine  COMPUTEEF  which  computes  the 
R  and  6  directed  components  of  the  E-field  (recall  that  there  is  no  (|)  component).  It  then 
calls  the  subroutine  CONVERTEF  to  transform  these  components  into  x,  y,  and  z  directed 
components  for  comparison  to  the  FDTD  results  that  are  in  Cartesian,  not  polar 
coordinates.  The  program  then  writes  out  the  value  of  the  E-field  to  a  data  file  for  Madab  to 
read. 

This  field  distribution  is  assumed  to  be  throughout  a  muscle-filled  space;  therefore, 
the  field  grows  everywhere  outside  the  sphere  of  muscle  tissue.  The  FDTD  code  defines 
the  E-Field  source  distribution  at  Rappaport's  maximum  9.45  cm  radius.  Defining  this 
surface  around  the  sphere  electrically  isolates  the  inside  and  outside.  In  this  manner,  the 
inside  of  the  sphere  should  contain  the  same  distribution  as  the  analytical  solution  within 
the  9.45  cm  radius.  Outside  this  radius  of  interest,  the  fields  will  be  drastically  different, 
thus  all  analysis  ignores  this  region. 

This  field  distribution  is  then  compared  with  the  FDTD  solutions  of  the  same  sphere 
geometry.  This  provides  an  analysis  to  test  the  validity  of  the  FDTD  code  and  its 
modifications.  This  simulation  provides  a  simple  duplication  of  Rapport’s  work. 

4.2  FDTD  Computer  Code 

The  code  for  producing  the  analytic  solution  is  used  to  specify  the  incident  field  in 
the  FDTD  space.  This  code  sets  up  an  electric  'shell'  around  the  sphere  by  setting  the  E- 
field  values  at  every  cube  along  the  9.45  cm  radius.  This  distribution  then  radiates,  and  the 
field  propagates  through  the  FDTD  grid  and  produces  the  same  results  as  the  analytic 
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version.  This  simulation  validated  the  use  of  FDTD  in  this  way  to  simulate  a  constant  E- 
field  around  a  sphere. 

4.2.1  FDTD  Geometry  Implementation 

Before  discussing  the  implementation  of  the  equations,  it  is  important  to  discuss  the 
definition  of  geometries  in  FDTD  and  the  layout  of  the  Yee  cell  [9].  An  understanding  of 
this  geometry  is  essential  to  understand  the  FDTD  equations. 

The  design  of  the  Yee  cell  offsets  the  magnetic  and  the  electric  fields  by  one  half  a 
cell  size.  Along  with  being  offset  by  a  half  step  in  space,  the  code  also  updates  the  E  and  H 
fields  a  half  time  step  apart.  Figure  4.2  shows  the  geometry  of  the  Yee  ceU.  The  cell 
shows  the  location  of  each  of  the  E  and  H  field  components.  The  user  creates  a  geometry 
by  assigning  a  material  ID  number  to  each  field  component  location  that  defines  its  electrical 
characteristics. 
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[9]. 


The  user  builds  FDTD  structure  by  assigning  a  material  ED  number  to  each  position. 
These  numbers  run  from  0  to  the  number  of  materials  in  the  space.  Zero  specifies  free 
space,  1  specifies  perfect  electric  conductors  (PEC)  and  numbers  two  and  up  have  a 
corresponding  permeability,  permittivity,  conductivity  and  magnetic  conductivity  stored  in 
a  matrix.  This  building  process  defines  six  matrices  containing  an  ID  number 
corresponding  the  material  at  each  position.  The  matrix  IDONE  contains  the  material  ID 
number  for  each  of  the  X  directed  components  in  the  FDTD  space;  likewise,  IDTWO  and 
IDTHRE  correspond  to  the  material  ID  number  at  the  Y  and  Z  locations  respectively. 
IDFOUR,  IDFIVE,  and  IDSIX  contain  the  ID  numbers  for  the  materials  located  at  the  X, 
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Y,  and  Z  components  of  the  magnetic  grid  in  the  FDTD  space.  These  six  matrixes  define 
the  objects  in  the  FDTD  space  [9]. 

4.2.2  FDTD  Implementation  in  FORTRAN 

Implementing  FDTD  in  a  computer  code  involves  solving  eq(15)  and  eq(16)  over 
time  and  space.  This  can  take  an  immense  amount  of  time  if  the  problem  space  is  large  or 
the  frequency  high.  Therefore,  it  is  important  to  implement  the  code  in  the  most  efficient 
manner  possible. 

The  first  step  in  doing  this  is  to  note  that  the  code  does  not  have  to  recalculate  the 
multipliers  before  each  field  component  for  every  cell  at  every  time  step  [9].  Since  most 


problems  will  have  a  small  number  of  different  scattering  materials,  the  program  calculates 
and  stores  these  multipliers  away  ahead  of  time  in  a  matrix.  To  do  this,  define  the 


following  FORTRAN  equalities  [9]; 
ECRLY(M)  = 


DT 


ECRLZ(M)  = 


ESCTC(M)  = 


EINCC(M)  = 


EDEVCN(M)  = 


(EPS(M)  +  SIGMA(M)  *  DT)  *  DY 

_ DT _ 

(EPS(M)  +  SIGMA(M)  *  DT)  *  DZ 
EPS(M) 


(EPS(M)  +  SIGMA(M)  *  DT) 
SIGMA(M)*DT 
(EPS(M)  +  SIGMA(M)  *  DT) 
DT*(EPS(M)-EPSO) 
(EPS(M)  +  SIGMA(M)  *  DT) 


At 

(e  +  (jAt)AY 
At 

(e  +  (TAr)AZ 
£ 

e  +  aAt 
aAt 

e  +  aAt 
{£-£jAt 
£+<7At 


Where  M  is  an  index  representing  the  different  material  types  present  in  the  problem  space. 
Defining  and  solving  these  equations  ahead  of  time  reduces  the  computation  time 
considerably.  These  numbers  are  then  input  into  eq  (19)  page  14  to  produce  a  form 
suitable  for  FORTRAN  code  [9]. 


EXS(I,J,K)  =  EXS(I,J,K)*  ESCTC(IDONE(I,J,K)) 

-EINCC(IDONE(I,  J,  K))  *  EXI(I,  J,  K) 

-EDEVCN(IDONE(I,  J,  K))  *  DEXI(I,  J,  K)  (6) 

+(HZS(I,  J,  K)  -  HZS(I,  J  - 1,  K))  *  ECRLY(IDONE(I,  J,  K)) 
-(HYS(I,  J,  K)  -  HYS(I,  J,  K  - 1))  *  ECRLZ(IDONE(I,  J,  K)) 


30 


Similar  equations  are  also  formulated  for  the  Y  and  Z  components  of  the  electric  field  as 
well  as  for  the  X,  Y,  and  Z  components  of  the  magnetic  field. 


Figure  4.3:  This  figure  shows  a  flow  diagram  of  the  FDTD  code. 


To  solve  for  the  field  over  time,  the  code  runs  through  the  time  steps,  recalculating 
these  equations  for  every  component.  For  each  of  these  loops,  the  code  first  runs  through 
all  the  X,  Y  and  Z  locations  of  the  electric  grid.  For  each  location  in  the  FDTD  space,  the 
X,  Y  and  Z- components  of  the  electric  field  are  updated  using  the  corresponding  equation 
similar  to  eq  (6).  The  program  then  calls  the  FEED  subroutine  that  forces  the  field  at  the 
feed  locations  to  be  equal  to  the  specified  source  wave  form.  The  time  is  then  incremented 
by  half  a  time  step  and  the  X,  Y  and  Z  components  of  the  magnetic  grid  are  updated  using 
the  corresponding  magnetic  equations  similar  to  eq  (6).  The  program  then  increments  time 
by  another  half  step  and  the  process  repeated.  Figure  4.3  shows  a  flow  diagram  of  this 
process. 

4.2.3  Code  Modifications 

Two  major  modifications  had  to  be  made  on  the  Penn  State  code.  First,  the  FEED 
subroutine  had  to  be  changed.  The  design  of  the  code  is  for  use  in  antenna  calculations; 
thus,  it  is  set  up  to  have  a  small  number  of  feeds,  one  for  each  antenna.  Forcing  an  E-field 
around  the  surface  requires  around  50,000  feeds.  Specifying  each  one  individually  proves 
to  be  both  time  consuming  and  memory  intensive.  To  avoid  this  pitfall,  the  FEED 
subroutine  is  modified  to  calculate  each  value  using  the  analytic  field  equations.  Figure  4.4 
and  figure  4.5  show  flow  diagrams  of  two  modifications  of  the  subroutine . 

The  first  modification  cycles  through  the  FDTD  space  and  determines  if  each  point 
is  along  the  surface  of  the  sphere.  When  it  comes  across  a  point  on  the  surface,  it  calls  the 
COMPUTEEF  subroutine  from  the  analytic  program.  This  computes  the  amplitude  and 
phase  for  the  R  and  THETA  components  of  the  E-field.  It  then  uses  a  simple  conversion 
routine  to  convert  into  x,  y  and  z  directed  field  components  and  forces  the  E-field  at  the 
location  to  have  the  desired  value.  The  Yee  cell  geometry  places  each  field  component  in  a 
slightly  different  location.  The  first  modification  computes  the  location  of  each  component 
and  then  figures  the  analytic  value  of  the  field  individually  for  each  component.  This 
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requires  computing  6  values  for  every  cube.  This  subroutine  repeats  every  time  step,  so 
increasing  its  speed  can  significantly  increase  computing  time. 

I  Sub.  Source.for  1 
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Figure  4.5:  Flowchart  of  FEED  subroutine  that  defines  an 
approximate  source  distribution  based  on  the  E-field  at  the  center  of 

each  cube. 
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To  reduce  the  computational  time  at  each  time  step,  the  second  modification 
approximates  the  source.  This  code  does  not  compute  a  different  value  for  each 
component.  Instead,  the  code  computes  the  analytical  value  the  components  would  have  at 
the  center  of  the  cube  and  gives  all  six  components  their  value  based  on  this  result.  A 
further  modification  to  decrease  computation  time  would  be  to  cycle  through  theta  and  phi 
at  the  radius  9.45  cm  rather  than  cycling  through  all  space  thereby  reducing  the  cycle 
length. 


The  second  major  modification  is  the  addition  of  a  peak  E-field  storage  array.  This 
allows  the  user  to  specify  an  area  over  which  the  peak  E-field  value  will  be  recorded. 
Figure  4.6  shows  the  code  of  the  implementation  of  this  in  the  DATSAV  subroutine. 

C  WRITES  OUT  MAX  VALUES  FOR  A  BLOCK  DEFINED  BY  TWO  POINTS 

C  (ONLY  USEFUL  FOR  SINGLE  FREQUENCY  RUNS) 

ELSE  IF  (NTYPE(NPT) .EQ.19)  THEN 
IF  (N.GT. ( .9*NSTOP) )  THEN 
REWIND ( 5 0+NPT) 

DO  246  I=IOBS(NPT) ,IOBSE(NPT) 

DO  247  J=JOBS (NPT) , JOBSE(NPT) 

DO  248  K=KOBS(NPT) ,KOBSE(NPT) 

EX=ABS(EXS(I, J,K)+FINC{I, J,K, 4) ) 

EY=ABS(EYS(I, J,K)+FINC(I, J,K,5) ) 

EZ=ABS(EZS(I, J,K)+FINC(I, J,K, 6) ) 

IF  (EX.GT.EXMd,  J,K)  )  EXM(I,  J,K)  =EX 
IF  (EY.GT.EYMd,  J,K)  )  EYMd,J,K)=EY 
IF  (EZ.GT.EZMd,  J,K)  )  EZMd,J,K)=EZ 
WRITE  (  (5 0+NPT)  ,  *)  EXMd,  J,K)  ,EYMd,  J,K)  ,EZMd,  J,K) 

248  CONTINUE 

247  CONTINUE 

246  CONTINUE 

ENDIF 


Figure  4.6:  The  Code  of  the  DATSAV  modification. 

This  is  only  useful  because  the  FDTD  simulations  are  using  one  electromagnetic  frequency. 
By  recording  the  peak  value  of  the  E-field,  the  power  deposited  can  be  easily  calculated 

IeP 

using  P  =  a.  Where  o  represents  the  conductivity  of  the  material.  Thus  a  simple 


conversion  into  the  frequency  domain  has  been  made  for  comparison  to  the  analytical 
results. 


35 


Simulation  description 

Cell  Size 

Srce  Radius 

Location  of  Schematic 

Homogeneous  Sphere 

2.55  mm 

9.45  cm 

N/A 

4-Layer  Laminated  Sphere 

2.55  mm 

9.45  cm 

Figure  4.6a 

3-Layer  Laminated  Sphere 

2.55  mm 

9.45  cm 

Figure  4.6b 

2-Layer  Laminated  Sphere 

2.55  mm 

9.45  cm 

Figure  4.6c 

Small  Head  w/  Muscle  Bolus 

2.55  mm 

9.45  cm 

Figure  4.7  and  4.8 

Small  Head  w/  Bone  Bolus 

2.55  mm 

9.45  cm 

Figure  4.7  and  4.8 

Large  Head  w/  Muscle  Bolus 

3.2  mm 

11.8  cm 

Figure  4.7  and  4.8 

Large  Head  w/  Bone  Bolus 

3.2  mm 

11.8  cm 

Figure  4.7  and  4.8 

Table  4.1:  This  table  gives  a  description  of  each  of  the  setups  along  with 

the  locations  of  their  schematics. 

These  two  modifications  allow  the  source  distribution  to  be  specified,  and  the 
resultant  field  to  be  recorded  and  analyzed.  The  first  modification  allows  the  FDTD  code  to 
implement  the  source  distribution  around  any  geometry  specified  in  the  FDTD  space.  The 
code  propagates  the  field  through  the  space  and  the  second  modification  allows  the 
frequency  response  to  be  recorded  for  analysis  in  matlab.  Appendix  B  gives  a  number  of 
matlab  routines  that  aid  in  analysis  of  the  field  distributions.  These  routines  make  it 
possible  to  look  at  any  slice  through  the  head  as  well  compute  errors  and  hot  spots. 

4.3  FDTD  RUNS 

Three  basic  forms  of  FDTD  simulations  are  done.  Table  4.1  gives  a  description  of 
each  setup  as  well  as  the  figure  numbers  of  the  schematics.  These  are  set  up  with  identical 
FDTD  space  parameters,  but  the  geometry  is  changed  for  each  simulation.  For  these 
simulations,  a  92x92x95  grid  is  set  up  with  cubic  cells  2.55  millimeters  on  a  side.  The 
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Courant  stability  condition  (page  16)  produces  a  49.19  picosecond  time  step.  The  outer 
sphere  in  each  simulation  has  a  radius  of  37  FDTD  cells,  or  9.45  cm.  The  analytical 
solution  sets  the  spherical  source  distribution  at  the  radius  of  37  FDTD  cells. 

The  first  simulation  is  on  a  simple  homogeneous  sphere.  This  simulation  serves  as 
a  validation  of  the  code.  It  is  to  demonstrate  the  ability  to  use  the  FDTD  code  with  an  entire 
incident  field  specified.  The  results  of  this  simulation  are  compared  with  the  analytic 
results. 

For  the  second  simulations,  a  series  of  laminated  spheres  simulates  the  skin-bone- 
brain  structure  of  the  human  head.  Figure  4.7  shows  a  schematic  of  each  of  these 
simulations.  The  first  of  these  is  a  two  level  sphere  the  outer  shell  being  bone,  the  inner 
being  muscle.  This  model  defines  an  inner  sphere  at  a  radius  of  32  FDTD  cells.  The  next 
simulation  is  of  a  three-layer  model  composed  of  an  outer  muscle  layer  to  simulate  the  skin 
and  bolus,  a  layer  of  bone  to  simulate  the  skull  and  an  inner  sphere  of  muscle  simulating 
the  brain.  The  radius  for  the  inner  sphere  was  18  cells,  the  outer  radius  for  the  bone  sphere 
is  19  cells  and  the  outer  sphere  is  again  at  a  radius  of  37  FDTD  cells.  A  four  layer  model 
comprises  the  final  simulation.  This  contained  an  outer  shell  of  bone,  to  simulate  a  bone¬ 
like  fluid  bolus,  a  thin  layer  of  muscle  tissue,  simulating  the  skin  layer,  a  bone  sphere  to 
simulate  the  skull  and  finally  an  inner  shell  of  muscle,  simulating  the  brain.  For  this 
simulation,  the  radius  of  the  inner  sphere  is  thirty-two  FDTD  cells.  The  outer  radius  of  the 
bone  layer  is  thirty-four  FDTD  cells,  and  the  radius  of  the  skin  shell  is  35  FDTD  cells. 
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34  cells 


Figure  4.7c:  2-Layer  Model 

Figure  4.7:  This  figure  shows  a  schematic  of  each  of  the  laminated 
sphere  simulations  (not  to  scale).  4.7a  is  the  four  layer  sphere, 
4.7b  is  the  three  layer  sphere  and  4.7c  is  the  two  layer  sphere. 
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An  FDTD  model  of  an  actual  human  head  provides  the  basis  for  the  final 
simulations.  This  model  is  made  by  converting  an  MRI  scan  of  a  human  head  into  a  head 
mesh  to  be  read  in  by  the  Penn  State  FDTD  code.  This  is  a  four  tissue  model  of  the  head. 
To  create  if,  the  MRI  scan  creates  a  map  of  the  different  materials  of  the  head.  It  is  then 
necessary  to  determine  the  actual  tissue  type  at  each  location.  Once  the  operator  determines 
the  exact  tissue  type,  the  typical  electrical  characteristics  of  each  tissue  are  assigned  to  each 
point  in  the  FDTD  space.  In  doing  this,  the  actual  permeabilities,  permitivities, 
conductivities  and  magnetic  conductivities  are  never  directly  measured.  Typical  values  are 
inserted  for  each  material  type.  This  provides  a  first  level  approximation  to  the  human 
head.  The  addition  of  a  basic  neck  extension  was  added  to  the  Penn  State  head  to  simulate 
the  inability  to  place  a  source  in  the  neck  region.  The  addition  of  a  neck  is  a  simple 
extension  of  the  final  layer  in  the  head  down  to  the  edge  of  the  FDTD  grid.  Figures  4.8 
and  4.9  show  two  cuts  through  the  head  model.  X=0  is  the  back  of  the  head,  Y=0  is  the 
left  side  of  the  head  and  7j=Q  corresponds  to  the  bottom  of  the  head. 
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X  Location 


Figure  4.8:  This  flgure  shows  a  horizontal  cut  through  the  center  of  the  head 
model.  X=0  corresponds  to  the  back  of  the  head  and  Y=0  corresponds  to  the 

left  side  of  the  head. 


A  liquid  bolus  surrounds  this  head  approximation.  This  bolus  allows  the  use  of  the 
spherical  source  distribution  on  the  edge  of  the  bolus,  providing  a  closer  match  in  material 
parameters  as  the  field  propagates  into  the  head.  Although  it  is  about  9.45  cm  to  the  center 
of  the  head,  it  is  not  spherical;  therefore,  the  full  size  head  cannot  fit  within  the  9.45  cm 
bolus  and  has  to  be  shrunk  for  the  initial  tests.  To  do  accomplish  this,  the  original  head 
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Coronal  Cut  Through  the  Head  Model  at  X=38.  Material  ID  # 
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Figure  4.9:  This  figure  shows  a  longitudinal  cut  through  the  center  of  the  head 
model.  Y=0  is  the  left  side  of  the  head  and  Z=0  is  the  bottom  of  the  head. 


mesh  cell  size  of  3.2  mm  is  reduced  to  2.554  millimeters.  The  material  characteristics  of 
the  bolus  can  be  set  at  any  value  desired.  The  first  simulation  gives  muscle-like 
characteristics  to  the  bolus.  The  second  simulation  uses  a  bone-like  bolus.  These  two 
simulations  provide  a  comparison  to  determine  which  material  characteristics  in  the  bolus 
are  preferred. 
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Along  with  the  simulations  using  a  small  head  model,  two  simulations  use  the  full 
size  human  head.  For  these  simulations,  the  cell  size  increases  from  2.554  millimeters  to 
3.2  millimeters  on  each  side  of  the  FDTD  cubes.  These  simulations  provide  an  opportunity 
to  determine  if  increasing  the  penetration  depth  in  the  actual  head  is  possible.  Brain  tissue 
is  less  lossy  than  muscle  tissue;  therefore,  it  is  conceivable  that  greater  penetration  may  be 
possible  in  the  presence  of  the  inhomogeneities  of  the  human  head.  These  two  simulations 
also  contain  one  simulation  with  a  bone-like  bolus  and  one  done  on  a  head  with  a  muscle¬ 
like  bolus. 
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5.  Results 


The  results  for  all  the  simulations  proved  to  be  very  encouraging.  The  analytic 
simulation  duplicates  Rappaport's  work  with  excellent  accuracy,  indicating  a  9.45  cm 
sphere  as  the  maximum  size  sphere  through  which  a  centrally  located  tumor  could  be 
irradiated.  The  FDTD  simulation  on  a  homogeneous  sphere  shows  excellent  correlation 
with  the  analytic  simulations.  These  results  show  it  is  possible  to  model  this  scenario  using 
the  FDTD  method.  The  laminated  sphere  simulations  demonstrate  a  "staircase"  effect  error 
introduced  from  the  cubed  approximation  of  spheres.  This  illustrates  an  error  to  watch  for 
in  the  head  simulations.  Finally,  the  head  simulations  prove  that  not  only  will  this 
technique  work  on  a  small  head  within  the  9.45  cm  optimal  radius,  it  wiU  also  work  on  a 
full  sized  human  head. 

All  the  figures  in  this  section  refer  to  levels  within  the  spheres.  These  levels 
correspond  to  grid  cells  within  the  FDTD  space;  e.g.,  X=38  refers  to  38th  plane  forward 
from  the  back  of  the  head.  FDTD  simulation  requires  a  10  cell  empty  space  region  around 
the  spheres  that  is  cut  out  to  clarify  the  results'  descriptions;  therefore,  the  levels  run  from 
1  to  75.  Since  the  spheres  were  37  cells  in  radius,  the  exact  center  of  the  spheres  and  head 
are  at  X=38,  Y=38,  Z=38.  The  bottom  of  the  head  is  at  Z=l,  the  back  of  the  head  is  at 
X=1  and  the  left  side  of  the  head  is  located  at  Y=l. 

5.1  Analytical  Simulation 

The  first  simulation  solves  the  analytic  solution  across  space  for  the  homogeneous 
sphere.  This  simulation  serves  as  the  control  to  ensure  FDTD  gives  reasonable  results  in 
this  treatment  scenario.  Figures  5.1  and  5.2  show  the  results  for  this  simulation.  These 
figures  show  a  cut  through  the  equator  of  the  sphere  (figure  5.1)  and  a  cut  through  a 
meridian  of  the  sphere  (figure  5.2).  Each  grid  increment  represents  an  increase  of  2.554 
millimeters.  Therefore,  the  37  cell  radius  sphere  is  equivalent  to  the  9.45  centimeter  radius 
source  shell.  The  matlab  analysis  removes  everything  outside  the  location  of  the 
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X  Location 


Analytical  Results  for  the  Homogeneous  Sphere  at  Z=38. 
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Figure  5.1:  This  Hgure  shows  the  analytical  results  for  the  power  distribution 
across  the  center  XY  plane  of  the  homogeneous  sphere  (normalized  to  one  at  the 
center).  Everything  outside  the  proposed  source  location  is  removed  for  clarity. 


theoretical  source  to  clarify  the  plots.  In  reality,  the  field  would  continue  to  grow  towards 
the  comers.  This  is  because  the  analytic  solution  describes  an  E-field  distribution  and  not  a 


radiating  source. 
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Z  Location 


Analytical  Results  for  the  Homogeneous  Sphere  at  X=38.  Normalized  Power 
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Figure  5.2:  This  figure  shows  the  analytical  results  for  the  power  distribution 
across  the  center  YZ  plane  of  the  homogeneous  sphere  (normalized  to  one  at  the 
center).  Everything  outside  the  proposed  source  location  is  removed  for  clarity. 

Figure  5.1  shows  the  central  XY  plane  of  the  sphere.  The  light  circle  in  the  center 
is  the  location  of  the  power  concentration.  Because  the  pattern  as  no  phi  dependence,  this 
pattern  is  completely  symmetric.  Figure  5.2  shows  the  central  YZ  plane  of  the  sphere. 

Again  the  light  spot  in  the  center  of  the  sphere  illustrates  the  power  focus.  Unlike  the  XY 
plane,  this  power  profile  is  not  symmetric.  This  would  be  expected  due  to  the  theta 
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dependence  produced  by  the  optimization  process.  Recall  that  the  final  field  had  very  little 
strength  near  the  poles  and  a  power  maximum  at  three  points  near  the  equator  (figure  3.6). 

5.2  FDTD  Simulations  on  a  Homogeneous  Sphere 

The  first  simulation  on  the  FDTD  sphere  proved  to  be  ineffective.  This  simulation 
produced  significant  spikes  along  the  poles.  An  analysis  of  this  simulation  displays  spikes 
of  over  four  hundred  times  the  value  of  the  E-field  at  the  center  along  the  top  and  bottom  of 
the  sphere.  Therefore,  the  power  is  not  completely  coupling  into  the  sphere.  Defining  the 
source  as  a  surface  causes  these  spikes.  That  is,  when  the  code  reaches  a  point  in  a  grid 
that  is  along  the  spherical  source,  it  sets  the  three  field  components  corresponding  to  that 
cell  location  to  the  value  the  optimal  analytical  solution  dictates.  This  produces  a  surface  of 
sources.  This  surface  contains  small  holes  along  the  stepped  edges  of  the  sphere.  Because 
of  these  holes,  the  source  does  not  isolate  the  inside  of  the  sphere  from  the  outside  and  the 
E-field  circulates  through  the  holes  which  prevents  it  from  completely  coupling  into  the 
muscle  sphere  and  creates  the  spikes.  This  obstructs  the  power  at  the  center  from  reaching 
its  full  potential  value. 

To  reduce  this  effect,  the  source  was  redefined  as  a  cubic  volume  instead  of  a 
surface.  In  other  words,  at  each  source  point  the  code  defines  six  source  components 
instead  of  three.  This  cubic  volume  eliminates  the  holes  and  creates  a  continuous  spherical 
source. 

There  remains  one  other  consideration  in  the  source  distribution.  Recall  that  the 
Yee  cell  geometry  places  each  field  component  in  a  slightly  different  spatial  location. 
Therefore,  to  be  completely  accurate  in  defining  the  source  distribution,  the  code  must 
recalculate  the  value  for  each  field  component  based  on  its  specific  position.  This  requires 
six  individual  calculations  for  each  source  cube  at  every  time  step,  adding  tremendous 
computation  time  to  an  already  slow  process. 

Therefore,  two  more  homogeneous  simulations  investigate  the  effect  of  an  exact 
source  distribution  versus  an  approximate  distribution  to  speed  computation  time.  The 
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approximate  source  code  calculates  what  the  source  component  values  should  be  at  the 
center  of  the  cube  and  defines  all  six  components  based  on  these  values.  This  reduces  the 
computation  time  by  a  factor  of  six  with  a  minimal  loss  in  accuracy.  Table  5. 1  shows  a 
comparison  of  the  results.  The  table  gives  the  percent  error  and  standard  deviation  versus 
the  analytic  results. 


FDTD  Simulation 

Mean  %  Error  Vs. 

Standard  Deviation 

Analytic  Simulation 

Vs.  Analytic 

Simulation 

Exact 

5.39% 

7.84% 

Approximate 

6.00%  . 

8.49% 

Tables. 1:  Comparison  of  exact  and  approximate  source 
distributions  for  homogeneous  sphere  FDTD  simulations 

The  error  shown  in  table  5.1  indicates  very  little  loss  in  accuracy  due  to  the 
approximation  of  the  source.  The  tremendous  gain  in  computation  time  and  small  loss  in 
accuracy  dictates  the  use  of  the  approximate  solution  for  remaining  simulations.  Figure  5.3 
shows  a  mesh  of  the  power  distribution  across  the  YZ  plane  of  the  FDTD  approximate 
source  simulation.  Again,  each  grid  increment  is  equivalent  to  2.554  millimeters.  The 
spikes  along  the  top  and  bottom  of  the  sphere  show  locations  where  the  field  is  not 
completely  coupling  into  the  muscle  tissue.  This  effect  is  pronounced  at  the  top  and  bottom 
of  the  sphere  because  of  the  polarization.  Since  the  tangential  component  of  the  E-field  is 
always  continuous  and  the  field  is  essentially  Z-polarized,  the  field  along  the  equator  is 
tangential  and  couples  completely.  The  field  near  the  poles  is  normal  to  the  sphere's 
surface.  Because  the  normal  E-field  is  discontinuous  by  a  factor  of  the  difference  in  the 
permitivities,  the  field  at  the  poles  does  not  completely  couple  into  the  muscle  tissue; 
forming  the  spikes. 
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E-Field  Strength  (V/m) 


E-Field  Across  the  Homnogeneous  Sphere  at  Z=38  (approximate  source). 
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Figure  5.3:  Mesh  of  the  power  distribution  across  the  YZ  plane  of  the 
homogeneous  sphere.  The  results  shown  are  of  the  FDTD  simulation  on  the 
homogeneous  sphere  with  an  approximate  source  distribution  (normalized  to  one  at 

the  center). 

Despite  this  incomplete  coupling,  the  field  profiles  show  good  correlation  to  the 
analytical  results.  Comparing  the  similarities  between  the  FDTD  field  solution  in  figure  5.4 
and  5.5  and  the  analytical  solutions  shown  in  figures  5.1  and  5.2  shows  this  correlation. 

Figure  5.4  is  a  cut  through  the  equator  of  the  sphere  showing  the  power  profile  calculated 
by  the  FDTD  model.  Figure  5.5  shows  the  power  profile  across  the  center  YZ  plane 
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through  the  sphere.  Both  of  these  plots  have  the  fields  set  to  zero  outside  the  sphere  to 
clarify  the  display.  Again  the  bright  spot  at  the  center  shows  the  focusing  of  the  power. 
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FDTD  Results  for  the  Homogeneous  Sphere  at  Z=38. 
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Figure  5.4 
sphere  with 


:  Power  profile 
an  approximate 


across  the  equatorial  cut  through  the  homogeneous 
source  distribution  (normalized  to  one  at  the  center, 


source  removed). 


Figures  5.6  and  5.7  shows  the  FDTD  power  profile  along  with  the  analytical 


solution  along  the  central  Y  and  Z  lines  in  figures  5.6  and  5.7.  These  plots  show  extremely 
strong  correlation  between  the  analytical  and  the  FDTD  solutions  along  the  central  lines. 
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Z  Location 


FDTD  Reults  for  the  Homogeneous  Sphere  at  X=38. 
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Figure  5.5:  Power  Proflle  across  the  central  YZ  cut  through  the  homogeneous 


sphere  with  an  approximate  source  (normalized  to 


one  at  the  center,  source 


removed). 


50 


Normalized  Power  (W/m'^2) 


Power  Comparison  Along  Line  X=38,  Z=38. 


—  FDTD  Approximation 

Figure  5.6:  Power  comparison  along  the  line  X=38,  Z=38  showing  Analytical 
and  FDTD  results  (normalized  to  one  at  the  center). 


51 


Normalized  Power  (W/m^2) 


Power  Comparison  Along  Line  X=38,  Y=38. 


Figure  5.7:  Power  comparison  along  the  line  X=38,  Y=38  showing  Analytical 
and  FDTD  results  (normalized  to  one  at  the  center). 
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Error  Between  Analytic  and  FDTD  Results  on  a  Homogeneous  Sphere  at  Z=38. 


Figure  5.8:  Equatorial  cut  showing  the  error  between  the  analytical  and  FDTD 

results  for  the  homogeneous  sphere. 


Figure  5.8  shows  the  error  between  the  analytical  simulation  and  the  FDTD 
simulation  with  an  approximate  source  distribution.  This  figure  shows  the  error  across  the 
equatorial  cut  of  the  sphere.  As  expected,  this  central  plane  has  exceptionally  low  error, 
with  a  maximum  around  8%.  With  such  strong  correlation  along  the  central  lines  and  the 
equatorial  cut  the  question  becomes  where  the  error  lies. 
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The  answer  to  this  question  comes  by  looking  at  figure  5.9  which  shows  the  error 
across  the  central  YZ  plane.  It  is  clear  from  this  mesh  that  the  majority  of  the  error  is 
coming  along  the  central  core  of  the  sphere  along  the  Z-axis.  This  is  expected  because 
figure  5.3  showed  earlier  that  the  power  was  not  completely  coupling  into  the  sphere  near 
the  poles.  Since  the  power  is  not  completely  entering  the  sphere  in  these  regions  due  to  the 
E-field  discontinuity,  the  power  profile  of  the  FDTD  solution  will  differ  from  that  of  the 
analytical  solution  near  the  poles.  In  addition,  the  interference  of  all  the  sources  in  the 
respective  hemisphere  forms  the  power  along  the  central  core.  Thus,  the  lack  of  complete 
coupling  at  the  poles  only  slightly  effects  the  power  near  the  edges.  However,  the  small 
effects  on  the  edges  will  combine  in  the  central  region  where  interference  occurs  producing 
expanded  error.  In  other  words,  although  the  majority  of  the  interference  occurs  at  the  ' 
focus,  there  will  also  be  some  amount  of  interference  in  the  central  region  that  will  combine 
the  small  loss  along  the  poles  into  a  significant  amount  of  error. 


Error  Between  Analytic  and  FDTD  Results  for  a  Homogeneous  Sphere  at  X=38. 


Y  Location  0  0  ^  Location 


Figure  5.9:  Central  YZ  cut  showing  the  error  between  the  analytical  and  FDTD 

results  for  the  homogeneous  sphere. 


5.3  FDTD  simulations  on  Laminated  Spheres 

The  next  series  of  simulations  are  done  on  laminated  spheres  to  investigate  the 
effects  of  adding  simple  inhomogeneous  geometries  into  FDTD.  These  simulations  all 
demonstrated  the  same  stair  step  spikes  introduced  by  the  stair  step  approximation 
necessary  in  modeling  a  sphere  as  a  series  of  cubes.  The  schematics  for  each  of  the 
simulations  is  shown  in  figrue  4.6. 
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E-Field  Strength  (V/m) 


Mesh  of  the  E-Field  Across  the  4-Layer  Sphere  at  Z=38. 


Figure  5.10;  3-D  plot  of  the  E-field  across  the  equatorial  cut  of  the  four-layer 

laminated  sphere. 

The  first  of  these  simulations  is  on  a  four-layer  sphere  (figure  4.7a).  The  material 
properties  of  the  outer  sphere  model  bone-like  material.  The  next  layer  down  patterns 
muscle  material.  The  third  layer  is  defined  as  bone  material.  Finally,  the  properties  of  the 
inner  sphere  imitate  muscle  tissue.  Therefore,  this  simulation  is  attempting  to  recreate  a 
liquid  bolus  with  bone-like  electrical  characteristics  around  a  head.  This  simulation  models 
the  head  as  a  muscle-like  brain  core  surrounded  by  a  spherical  'skull'  of  bone  and  finally  a 
thin  'skin'  layer  of  muscle  tissue. 


80 
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Figure  5.10  shows  the  E-field  distribution  across  the  central  XY  cut  of  the 
laminated  sphere  normalized  to  one  at  the  center.  The  field  along  the  equator  of  the  sphere 
is  predominately  directed  in  the  Z  direction.  Thus,  along  this  equatorial  cut,  the  field  will 
be  tangent  to  the  surface  of  each  sphere.  Because  the  tangential  E-field  has  to  be 
continuous  across  a  dielectric  boundary,  along  this  cut  the  field  should  be  continuous 
everywhere.  Figure  5. 10  shows  this  is  not  the  case.  There  are  E-field  spikes  along  the 
edges  of  the  inner  spheres.  The  stair  stepped  edges  of  the  spheres  causes  this 
phenomenon.  Instead  of  having  a  smooth  surface,  small  cubes  compose  the  edge.  At 
some  point  a  given  layer  of  cubes  ends  producing  an  edge.  This  edge  produces  a  locally 
horizontal  rather  than  vertical  surface.  At  this  point  the  E-field  becomes  normal  rather  than 
tangential.  The  amount  of  the  discontinuity  of  the  normal  E-fields  is  proportional  to  the 
difference  in  relative  permitivities  of  the  two  materials;  therefore,  the  field  spikes  in  the 
cube. 
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E-Field  Across  the  4-Layer  Sphere  at  Z=39.  E-Field  (V/m) 
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Figure  5.11:  2-D  image  of  the  E-field  distribution  across  the  Z=39  cut  of  the  4- 

layer  laminated  sphere. 


Figures  5.11  and  5.12  illustrate  the  effect  of  the  staircase  errors.  Figure  5.11 
shows  the  E-field  distribution  across  the  sphere  in  a  two-dimensional  image.  The  spikes 
show  up  as  single  pixels  slightly  darker  or  lighter  than  those  around  them.  Figure  5. 12  is 
plot  of  the  changes  in  the  IDTHRE  components  from  level  Z=39  to  level  Z=40.  Recall  that 
the  IDTHRE  array  contains  the  material  ID's  for  the  Z-directed  geometry  components. 
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Changes  in  IDTHRE  Components  of  4-Layer  Sphere  From  Z=40  to  Z=39. 
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Figure  5.12:  This  figure  shows  a  2-D  image  of  the  changes  in  the  Z-directed 


material  components  from  Z=39  to  Z=40. 


That  is,  where  the  IDTHRE,  or  Z-directed  component  changes  from  one  level  to  the  next, 
there  is  a  locally  horizontal  surface  where  it  should  be  vertical.  Therefore,  figure  5.12 
shows  only  the  location  of  each  stair  step.  Investigation  of  figures  5.11  and  5. 12  reveals 
that  every  spike  corresponds  to  the  edge  of  a  stair  step.  These  spikes  are  numerical  artifacts 
which  will  not  occur  in  actual  practice. 
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E-Field  Across  4-Layer  Sphere  at  X=38. 
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Figure  5.13:  This  figure  shows  a  3-D  mesh  of  the  E-Held  distribution  across  the 

central  YZ  cut  of  the  four-layer  model. 

Figure  5.13  shows  the  E-field  distribution  across  the  central  YZ  plane  of  the  4-layer 
sphere.  As  expected,  the  field  is  continuous  along  the  central  Y  line.  The  discontinuities 
through  the  center  spheres  become  more  apparent  near  the  poles.  This  is  because  the  field 
is  predominately  Z-directed;  therefore  nearer  the  poles  the  field  is  normal  rather  than 
tangential  to  the  dielectric  interface  so  the  discontinuity  will  be  larger. 


60 


E-Field  (V/m) 


The  next  simulation  incorporated  three  layers  as  compared  to  the  former  4-layer 
model.  The  3-layer  simulation  produced  the  same  stair  step  spikes  along  with  a  new 
phenomenon .  The  3-layer  model  had  a  much  smaller  'skull'  layer  than  any  of  the  other 
models  (figure  4.7b).  As  the  field  nears  the  center,  the  radial  component  becomes  more 
pronounced.  Thus,  along  the  equatorial  cut,  shown  in  figure  5.14,  the  field  has  a 
noticeable  discrepancy  as  it  passes  through  the  dielectric  interfaces.  The  radial  component 
is  more  prominent  near  the  center  of  the  sphere  causing  the  discontinuity.  The  field  thus 

has  a  larger  normal  component  to  produce  the  discontinuity. 

E-Field  Across  2-Layer  Sphere  at  Z=38. 


Figure  5.14:  This  Hgure  shows  the  E-Field  across  the  equatorial  cut  of  the  3-layer 
sphere  model  (normalized  to  one  at  the  center). 
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E  field  along  line  X=38,  Z=38. 
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Figure  5.15:  Plot  of  the  X  and  Y  E-fleld  components  along  an  equatorial  radial 
line  in  the  3-layer  Sphere  model  (E-field  normalized  to  one  at  the  center). 

Figure  5.15  is  a  plot  of  the  x  and  y  components  of  the  E-field  along  an  equatorial 
radial  line.  The  Z  directed  component  contains  the  theta  component  along  equatorial  plane. 
The  X  and  Y  directed  components  contain  the  radial  component.  Figure  5.15  shows  these 
two  components  along  an  equatorial  line.  From  this  plot  the  discontinuity  can  be  seen  to  be 
due  to  the  larger  radial  component  that  is  normal  to  the  interface. 

Using  the  2-layer  model  (figure  4.7c)  in  the  next  simulation  showed  the  same  stair 
step  spikes  seen  in  the  other  laminated  sphere  simulations.  Figures  5.16  and  5.17  show 
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the  results  of  this  simulation.  Figure  5.16  is  a  mesh  of  the  field  strength  across  the 
equatorial  cut  of  the  2-layer  sphere.  Figure  5.17  shows  the  field  strength  across  the  central 
YZ  plane. 

E-Field  Across  2-Layer  Sphere  at  Z=38. 


X  Location 


0  0 


Y  Location 


Figure  5.16:  This  Hgure  shows  a  3-D  mesh  of  the  E-field  strength  across  the 
equatorial  plane  of  the  two-layer  model  (normalized  to  one  at  the  center). 
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E-Field  Across  2-Layer  Sphere  at  X=38. 


Figure  5.17:  This  figure  shows  a  3-D  mesh  of  the  E-Held  strength  across  the 
central  YZ  plane  of  the  two-layer  model  (normalized  to  one  at  the  center). 

These  simulations  on  laminated  spheres  show  several  important  aspects  that  will  be 
seen  in  the  simulations  on  the  actual  head  model.  The  most  important  of  these  lessons  is 
the  presence  of  the  "staircase"  spikes.  The  head  simulation  shows  these  spikes  as  well.  It 
is  important  to  keep  in  mind  that  extended  spikes  due  to  a  true  horizontal  surface  will  occur; 
the  single  cell  spikes,  however,  are  due  to  the  "staircase"  error  and  would  not  occur  in 
actual  practice.  Another  lesson  this  illustrates  is  the  spiking  due  to  the  horizontal  surfaces 


64 


near  the  poles.  Computing  the  power  will  significantly  reduce  the  spiking  due  to  the  very 
low  conductivities  of  the  bone  tissue  where  the  spiking  occurs.  If  it  does  become  a 
problem,  going  back  to  the  original  field  distribution  that  concentrates  the  source  power 
more  near  the  equator  can  reduce  this  spiking. 

5.4  FDTD  Simulations  on  the  Head  Model 

After  investigating  the  effects  of  the  simple  inhomogeneities,  The  next  simulation 
places  the  head  model  derived  from  an  MRI  scan  within  the  homogeneous  sphere.  The 
sphere  will  then  act  as  a  model  of  a  liquid  bolus  placed  around  the  patient's  head.  This 
bolus  is  used  to  help  the  field  couple  more  completely  into  the  tissue.  To  investigate  the 
effects  of  bolus  composition,  two  different  bolus  types  are  used,  a  bone-like  bolus  and  a 
muscle-like  bolus.  This  provides  the  basis  for  five  different  head  simulations.  The  first  is 
an  error,  reported  here  to  demonstrate  the  importance  of  proper  characteristics.  The  next 
two  simulations  use  a  shrunken  head  with  both  a  bone-like  and  a  muscle-like  bolus.  The 
final  two  simulations  use  a  full  size  head  model,  one  with  a  bone-like  bolus  and  one  with  a 
muscle-like  bolus. 

The  first  simulations  on  the  head  model  produce  troubling  results.  The 
inhomogeneous  head  model  had  absolutely  no  correlation  to  the  simulation  on  the  head 
model.  A  mesh  across  the  central  horizontal  slice  is  shown  in  figure  5.18.  This  plot 
shows  tremendous  spiking  throughout  the  edge  regions  and  practically  no  focusing. 
Further  investigation  shows  these  discrepancies  are  due  to  inaccurate  material 
characteristics  in  the  head  model.  The  model  had  been  given  the  material  characteristics  of 
fat  rather  than  muscle  for  the  brain  tissue.  These  results  demonstrate  the  importance  of  the 
optimization  process.  The  focusing  could  be  accomplished  on  such  a  head;  however,  a 
new  optimization  would  have  to  be  done  to  account  for  the  different  material  characteristics 
associated  with  fat. 


65 


Power  Across  Fat  Head  Model  at  Z=38. 


Figure  5.18:  Power  across  the  central  horizontal  slice  of  the  head  model 
(normalized  to  one  at  the  center).  These  results  are  of  a  head  with  a  bone*like 

bolus  and  brain  tissue  defined  as  fat. 


The  remaining  four  simulations  again  use  the  model  of  the  head  developed  from  an 
MRI  scan,  except  with  proper  material  characteristics.  The  material  characteristics  define 
three  different  materials  with  muscle-like  properties  representing  skin  and  white  and  gray 
matter,  and  one  material  with  the  low-water  content  properties  of  bone.  The  first  two  of 
these  uses  a  head  that  is  reduced  to  fit  within  the  9.45  cm  source  shell.  The  second  two  use 
a  full  size  head,  with  an  enlarged  source  shell  to  fit  around  the  enlarged  head.  Each  of 
these  sets  contains  one  simulation  with  a  bone-like  bolus  and  one  done  with  a  muscle-like 
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Power  {W/m''2) 


bolus.  Figures  5.19  through  5.22  show  a  plot  across  the  central  XY  plane  of  each  of  these 
simulations. 

Horizontal  cut  of  the  small  head  with  muscle-like  bolus  at  Z=38. 


Figure  5.19:  Power  across  the  central  XY  plane  of  the  simulation  of  a  small  head 
model  with  a  muscle  bolus  (normalized  to  one  at  the  center). 
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Horizontal  cut  of  the  large  head  with  muscle-like  bolus  at  Z=38. 


X  Location 


Y  Location 


Figure  5.21:  Power  across  the  central  XY  plane  of  the  simulation  of  a  large  head 
model  with  a  muscle  holus  (normalized  to  one  at  the  center). 
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Power  (W/m^2) 


Horizontal  cut  of  the  large  head  with  bone-like  bolus  at  Z=38. 


Figure  5.22:  Power  across  the  central  XY  plane  of  the  simulation  of  a  large  head 
model  with  a  bone  bolus  (normalized  to  one  at  the  center). 


Table  5.2  shows  a  comparison  of  the  analysis  statistics  on  each  of  the  simulations. 
Center  mean  is  a  measure  of  the  mean  power  in  the  center  three-cell  cube  of  the  head. 
Because  the  power  is  normalized  to  one  at  the  center,  this  will  always  be  near  one.  The 
smaller  it  is,  however,  the  more  narrow  the  spike  will  be  near  the  center,  thereby  producing 
a  more  accurate  focusing.  Total  mean  is  a  measure  of  the  mean  power  everywhere  else  in 
the  head.  This  does  not  include  the  power  in  the  bolus  because  circulation  of  the  bolus 
liquid  can  prevent  overheating  in  the  bolus.  The  lower  this  value,  the  smaller  the  risk  will 


70 


be  of  over-heating  undesirable  locations.  The  third  value  in  the  table  is  the  ratio  of  the 
mean  power  at  the  center  to  the  total  mean  power.  The  higher  this  value,  the  easier  it  wiU 
be  to  heat  the  tumor  without  affecting  healthy  tissue.  The  final  three  numbers  are 
percentages  of  cells  over  various  thresholds,  the  first  column  gives  the  percent  of  cells 
with  a  value  greater  than  1.  Each  of  these  points  will  reach  temperatures  higher  than  the 
central  tumor,  thus  destroying  the  healthy  tissue  at  that  location.  It  is  preferable  this 
number  remains  small.  The  next  two  are  similar  to  this,  however  they  are  the  percent  of 
cells  over  .9  and  .8  the  power  at  the  center.  Although  these  points  would  not  be  raised  to 
fatal  temperatures,  the  high  temperature  levels  may  be  dangerous  to  the  healthy  tissue. 


Simulation 

Cntr 

Total 

CM/TM 

%  cells 

%  cells 

%  cells 

Mean 

Mean 

over  1 

over  .9 

over  .8 

Fat  Head  with  Bone 

.9910 

2.946 

.3364 

72.7  % 

75.9% 

79.3% 

Bolus 

(W/m^2) 

(W/m''2) 

Small  Head  with 

.8982 

.2187 

4.108 

.26% 

.39% 

.58% 

Bone  Bolus 

(W/m^2) 

(W/m^2) 

Small  Head  with 

.8795 

.1142 

7.203 

.01% 

.05% 

.08% 

Muscle  Bolus 

(W/m^2) 

(W/m''2) 

Full  Head  with  Bone 

.9057 

.2942 

3.078 

1.66% 

2.33% 

3.32% 

Bolus 

(W/m'^2) 

(W/m'^2) 

Full  Head  with 

.8436 

.2389 

3.531 

1.33% 

2.0% 

2.89% 

Muscle  Bolus 

(W/m^2) 

(W/m'^2) 

Table  5.2:  Statistics  on  the  FDTD  Head  model  simulations 

Table  5.2  shows  the  muscle  bolus  is  the  preferred  bolus  type.  However,  as  the 
head  size  becomes  larger,  the  bolus  composition  becomes  less  important.  This  is  because 
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the  increased  power  necessary  in  the  larger  head  models.  The  spikes  are  predominately 
introduced  in  the  neck  region.  This  is  because  the  source  distribution  is  abruptly  turned  off 
at  the  edge  of  the  neck.  Because  there  is  an  abrupt  change  in  the  source  distribution  from 
outside  and' inside  the  neck  region,  the  field  diffracts  along  the  edge  of  the  source 
distribution  introducing  the  large  spikes.  The  larger  this  discrepancy  the  larger  the 
diffracted  field  spikes.  In  the  small  head  models,  the  change  in  material  characteristics  in 
going  from  the  bone  of  the  bolus  to  the  muscle  of  the  neck  has  the  predominate  effect  in 
introducing  the  diffracted  field.  In  the  larger  head,  the  higher  source  power  makes  the 
stepped  edge  of  the  source  power  the  larger  contributor  to  the  diffracted  field.  Figure  5.19 
shows  the  power  spikes  in  the  neck  for  each  of  the  four  simulations. 

In  both  cases,  the  muscle  bolus  produces  a  more  precise  heating  pattern.  This  is  • 
because  the  optimization  is  designed  around  a  muscle  sphere.  The  phase  of  the  source  is 
designed  for  propagation  through  muscle.  Passing  through  bone  changes  the  velocity  of 
the  wave,  thereby  reducing  the  focusing,  and  smearing  the  central  spike  across  a  larger 
area. 
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6.  Conclusions 


This  research  provided  seeral  valueable  insights  into  electromagnetic  therapy  as  well 
presenting  some  important  areas  for  further  research.  The  final  simulations  on  the  complete 
head  model  illustrated  the  potential  of  microwave  therapy.  These  simulations  showed  with 
proper  optimization,  it  may  be  possible  to  irradiate  a  deep  set  tumor  within  a  full  size 
human  head  with  adequate  precision.  The  research  also  demonstrated  the  need  for 
improved  head  modeling  to  ensure  proper  material  characteristics  are  given  to  the  tissue  of 
the  head. 

6.1  Findings 

FDTD  proved  to  be  an  effective  tool  in  this  field  of  research.  The  initial  simulations 
duplicating  the  optimal  analytical  solution  derived  by  Rappaport  demonstrated  the  validity 
of  the  FDTD  numerical  approach  for  treatment  scenario. 

The  inhomogeneous  simulations  on  laminated  spheres  demonstrated  the  effects 
inherent  in  an  FDTD  Approach.  The  stairstep  approximation  of  a  spherical  surface 
introduced  'stairstep  spikes'  due  to  a  locally  horizontal  surface  along  an  otherwise  vertical 
interface.  In  addition  to  these  spikes,  the  FDTD  approach  proved  to  be  difficult  in 
obtaining  complete  coupling  from  the  source  to  the  sphere  of  treatment.  These  drawbacks 
not  only  demonstrate  the  effects  of  an  FDTD  approximations,  but  also  introduce  an  area  to 
be  considered  in  microwave  treatment.  Because  the  head  is  not  spherical,  there  will  be 
horizontal  sections  of  interface  along  what  would  otherwise  be  a  vertical  interface.  For 
example,  along  the  base  of  the  mandible,  there  is  a  long  vertical  interface.  In  addition, 
along  the  inside  of  the  occipital  cavity  there  are  vertical  interfaces.  All  of  these  areas  wiU  be 
candidates  for  the  'stairstep  spiking'  phenomenom. 

The  simulations  on  the  actual  head  models  proved  very  encouraging.  The  results 
demonstrated  a  muscle-like  bolus  is  preferred  over  a  bone-like  bolus.  The  muscle-like 
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bolus  provided  better  precision  as  weU  as  reduced  spiking  outside  of  the  center.  The 
source  distribution  placed  around  a  small  head  within  a  spherical  water  bolus  produced 
excellent  focusing.  There  were  only  20  cells  above  that  at  the  center,  aU  of  which  were 
confined  to  the  neck  region.  This  came  about  due  to  the  stepped-edge  of  the  source 
distribution  near  the  neck.  There  was  no  taper  down  as  the  source  neared  the  edge  of  the 
neck,  thus  a  diffracted  field  was  created  near  the  edge  of  the  source  distribution  creating 
several  small,  but  strong  spikes  in  the  neck  region. 

These  spikes  are  troubling,  although  not  discouraging.  As  mentioned  above,  there 
was  no  attempt  made  to  taper  the  edge  of  the  field  distribution,  in  addition  the  neck  region 
itself  was  extremely  approximate.  The  MRI  scan  produced  no  neck  region,  therefore  the 
head  was  simply  continued  down  from  the  base  to  the  edge  of  the  FDTD  space.  There  was 
no  spinal  chord,  esophagus,  or  trachea  modeled.  In  addition,  the  head  itself  was 
approximated  from  an  MRI  scan.  Although  serving  to  produce  an  example  of  what  could 
be  possible  in  this  area,  it  still  needs  a  great  deal  of  research  and  refinement. 

The  final  head  simulations  showed  it  is  possible  to  get  good  resolution  at  the  center 
of  an  actual  human  head  with  minimal  destruction  to  healthy  tissue.  With  the  addition  of  an 
optimization  process,  this  precision  can  be  furhter  enhanced.  This  can  also  be  used  to 
reduce  the  unwanted  fatal  spikes  in  the  neck  region. 

The  simulations  done  on  a  full  size  head  demonstrated  two  critical  lessons.  First, 
these  are  the  first  simulations  to  show  FDTD  can  be  used  on  a  complete  3-D  model  of  a 
human  head  in  a  treatment  scenario.  Second,  and  most  important,  these  simulations  show 
that  even  without  any  form  of  optimization  to  account  for  the  inhomogeneous  structure  of 
the  human  head,  it  is  possible  to  irradiate  a  deep-set  tumor  with  decent  precision. 

6.2  Recommendations 

The  results  obtained  in  this  research  shows  an  avenue  for  continued  research  in  the 
field  of  microwave  therapy.  The  simulations  on  the  head  model  demonstrated  the  ability  to 
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concentrate  microwave  energy  in  an  actual  head  model.  The  results  obtained  thus  far 
contain  no  optimization  done  based  on  the  inhomogeneities  of  the  head. 

Several  areas  present  themselves  as  prime  candidates  as  a  result  of  the  simultations 
presented  here.  The  most  intuitively  obvious  of  these  follow-on  oppurtunities  is  to  derive  a 
optimal  source  distribution  taking  into  account  the  presence  of  the  neck.  The  inability  to 
place  the  source  within  the  neck  region  requires  some  form  of  optimization  in  defining  the 
taper.  This  can  be  done  in  several  ways,  by  introducing  a  spatial  filter  or  by  customizing 
the  spherical  harmonics  to  concentrate  the  power  near  the  top  reducing  it  near  the  lower 
pole. 

A  comparison  of  the  FDTD  results  to  actual  measurements  of  a  real  head  also  needs 
to  be  done.  This  research  demonstrated  FDTD's  ability  to  duplicate  analytic  results.  It  was 
assumed  the  ability  to  model  analytic  homogeneous  solutions  implied  the  ability  to  model 
inhomogeneous  geometries.  A  simulation  of  these  solutions  compared  to  actual 
measurements  is  needed  to  further  prove  this  method  is  a  valid  modeling  tool. 

Another  area  open  for  research  is  in  developing  a  more  accurate  modeling  system 
for  creating  the  head  mesh.  The  head  mesh  used  in  this  research  is  only  an  approximate 
rendering  of  an  MRI  scan.  The  method  used  to  model  the  brain  has  no  connection  to  the 
actual  electrical  charateristics  of  the  that  individual's  head  tissue.  The  parameters  used  were 
'typical'  values  for  average  tissue.  The  actual  properties  in  real  tissue  vary  immensely  from 
individual  to  individual.  For  this  reason,  a  system  is  needed  to  non-invasively  measure 
material  parameters  of  the  head  region. 

This  research  serves  to  illustrate  the  potential  inherent  in  microwave  treatment.  The 
results  show  that  electromagnetic  field  propagation  can  be  predicted  in  an  inhomogeneous, 
three-dimensional  scenario.  In  addition,  the  resutls  show  that  good  penetration  and 
precision  can  be  obtained  in  an  actual  human  head  with  proper  optimization.  This  research 
clearly  indicates  the  need  for  continued  research  in  the  area  of  microwave  hyperhtermia  of 
cancerous  tissue. 
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Appendix  A:  The  Vector  Wave  Equation  in  Spherical  Coordinates 

The  source  optimization  developed  by  Rappaport  and  Morganthaler  is  based  upon 
solution  to  Ae  wave  equation  in  the  spherical  coordinate  system.  The  geometry  of  this 
coordinate  system  is  shown  in  figure  2.1.  Before  discussing  this  optimization  it  is 
important  to  go  over  the  basics  of  the  solution  to  the  spherical  vector  wave  equation. 


Figure  A.l;  The  Spherical  coordinate  system  as  used  in  this  thesis 

Solution  to  the  Spherical  Vector  Wave  Equation 

Any  electromagnetic  wave  in  an  isotropic,  homogeneous,  source  free  region  must 
satisfy  the  vector  wave  equation  [10] 
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Where  C  represents  a  field  vector  (E,  B,  D,  or  H)  and  |i,e,  and  a  represent  permeability, 
permittivity  and  electrical  conductivity  as  usual.  The  analysis  presented  here  will  deal 
exclusively  with  the  electric  and  magnetic  field  intensities,  E  and  H.  The  linearity  of  this 
equation  allows  the  vectors  E  and  H  to  be  expanded  by  Fourier  analysis  into  time 
harmonics  of  the  form  C  =  with  no  loss  in  generality.  This  allows  the  entire  time 

dependence  to  be  expressed  as  and  will  be  left  out  in  the  interest  of  clarity.  If  it  is 

then  noted  that  V^C  =  VV  •  C  -  V  x  V  x  C  and  that  eiico^  +  iaiico  =  ,  eq  (1)  can  be 

rewritten  as 

VV-C-VxVxC  +  A:^C  =  0.  (2) 

There  exists  a  solution  set  to  equation  (2)  containing  three  forms  of  waves,  a 
longitudinally  polarized  wave  and  two  transverse  waves  commonly  referred  to  as 
transverse  electric  (TE)  and  transverse  magnetic  (TM)  waves  dependent  upon  whether  the 
magnetic  or  electric  field  is  tangential  to  the  direction  of  propagation  (in  this  case  along  a 
radial  line  from  the  center  of  the  sphere  of  interest).  This  thesis  is  based  only  on  the  electric 
field  solution  for  a  TM  wave. 

To  find  an  expression  for  the  TM  wave,  we  will  look  for  a  magnetic  field  solution 
to  (2)  of  the  form 

H  =  Vx(fM(/?)v^)  (3) 

Where  m(R)  is  a  function  of  R  yet  to  be  determined  and  V'^is  a  scalar  function  of  R,  0,  and 
0  and  a  solution  to  the  equation 

VV  +  ^V  =  0-  (4) 


Employing  the  del  operator  of  the  spherical  coordinate  system,  the  three  components  of  H 
can  be  written  as 


H,=0, 


H«  = 


Rsmdd(l>^  ^  RdG^ 


(5) 


The  divergence  of  the  cross  product  of  any  two  vectors  is  identically  zero,  thus  the 
divergence  of  H  will  be  zero  and  (2)  can  be  rewritten  as 

VxVxH-/t'H  =  0.  (6) 

Expanding  this  out  in  spherical  coordinates  reveals  that  the  f  component  is  identically  zero 


and  the  6  component  is  satisfied  by  the  equation 


n  sin^  d  I  ddd0  dd 


sind—iuii/)  + 


1 

Rhm^ed(l)^ 


iuV)  + 


RsinO  dr  d(j> 


f  .  d  ,  , 


R  de 


The  0  component  is  satisfied  by  the  equation 

d  1  \  d  .  ^  d  .  ,1  d  1 

ddR^smdldS  dd  ^  J  dOR^&ix 


J _ ^ 

sin^  6  d(^ 


Rdr^dO' 

Both  eq(7)  and  eq(8)  will  be  satisfied  by  the  equation 


1  d^  ,  ^  d  ,  , 

t3j3s<“V')+-:-3s("V') 


R  de 


{UW)  +  —:: - -  Smd  —  (UW)  +— r - ,  -i  . 

R^ sin 6  d6\_  dd^  ]  sin^  d  d(j)‘ 


■{uy/)  +  k^u\i/  =  0.  (9) 


If  we  know  take  a  look  at  the  requirement  placed  on  y/by  eq  (4);  expanding  that  equation 


yields 


1  df  2  1  df  .  1 

— r—  r  — +—5 - -  sm0— ^  +— 5 — 5 ^ 

R^  dry  dr  J  R^  sin  0  ^0  V  dd  J  R^  sin^  d  d<l>^ 


+^V=o. 


When  this  is  done,  it  immediately  becomes  apparent  that  if  we  let  m(R)=R  eq.  (5)  will  be 


satisfied. 


This  shows  one  solution  to  the  vector  wave  equation  is  given  by  H  whose 


components  are 


H,  =  0,  He=— 

sin0  d<t> 


H.  = 


dyf 


The  optimization  is  based  upon  the  electric  field  solution  to  a  TM  wave  which  can  easily  be 


found  from  the  relation 


jtE  =  VxH.  (12) 


Which  produces  an  electric  field  of  the  form 
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(13) 


1  d\R\ir) 
k  dR^ 


+  k^R\ir, 


1  d^{R\(f) 
Tr  dRdd  ’ 


1  d\R\l/) 

kRsmd  dRdp 


Solution  for  the  Elementary  Wave  Function 


Those  paying  attention  will  be  quick  to  point  out  that  this  expression  for  the  electric 
field  means  absolutely  nothing  without  an  expression  for  yf.  You  will  remember  this 
variable  was  introduced  as  the  solution  to  the  scalar  wave  equation  (4)  in  spherical 
coordinates.  In  order  to  find  this  solution,  start  with  the  assumption  that  it  can  be  written 
as 

xi/  =  f(R,d,(t>)e-“^.  (14) 


Which  will  result  in  no  loss  of  generality  if  the  time  dependence  is  left  out  as  described 
earUer. 


Therefore,  in  a  homogeneous  isotropic  medium,  the  function/ must  satisfy  the 

scalar  wave  eq  (4)  which,  when  expanded,  looks  like  eq  (10).  If/can  be  written  as 
/  =  /fi(^)/0(^)/« (0)  >  then  eq(lO)  is  separable  into 

dR 

AY 


R"  +  2R^  +  {k^R^  -  )/«  =  0 


1 


dR^ 


,2  A 


sin  6  dd\ 


sin0 


ddj 


sin  6 


+^V,=o 


d<t>^ 


(15) 

\fe=0  (16) 

(17) 


Where  p  and  q  are  introduced  as  variables  of  separation  yet  to  be  determined.  Since  we  are 

working  in  a  homogeneous,  isotropic  medium,  the  properties  of  the  medium  will  be 

independent  of  equatorial  angle  (j) .  In  addition,  it  is  physically  impossible  for  a  field  to 

have  two  different  values  at  the  same  point  in  space,  therefore  eq(17)  must  be  periodic 

function  of  (])  with  a  period  of  27C.  To  avoid  the  complications  involved  in  introducing  a 

cos  m<t> 

complex  function,  set  f.  =  .  Thus,  q  must  be  an  integer  of  the  form 

^  sinm0 


m=0,±l,±2,.... 
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Finding  a  value  for  q  is  slightly  more  complicated.  To  find  this  solution,  start  by 

noting  the  remarkable  resemblance  of  eq  (16)  to  the  associated  Legendre  equation  [11] 

+  +  =  O  (18) 

Sind  day  dd  J  \_  sin  0_ 

In  fact,  if  we  set  =  n(n  + 1),  where  n=0,l,2...  we  have  the  same  equation,  if  fg  is  an 

associated  Legendre  polynomial  of  the  form 


fe(^)  =  P^^icosd)  =  (l  -  cos^  dj 


dcos"'  d 


P„(cosd).  (19) 


This  all  works  out  mathematically  but  it  is  also  important  to  understand  it 
physically,  fg  must  be  periodic  in  theta  with  a  period  of  2%,  symmetric  about  K,  single 


valued  and  finite  for  this  to  physically  make  sense.  Since  we  are  defining  the  Legendre 
function  as  a  function  of  cosd,  it  will  be  periodic  in  theta  with  a  period  of  27C  and 


symmetric  about  tc.  Examining  the  definition  in  eq(19)  it  quickly  becomes  apparent  that  the 
associated  Legendre  function  will  be  zero  at  ±1,  and  single- valued  and  finite  everywhere  in 
between.  With  this,  p  and  q  have  been  determined  as  well  as  solutions  for 

fg(d)  mdf,((t>). 

The  final  step  is  to  determine  a  solution  for  fifR) .  This  stems  from  the  solution 
to  eq  (15)  with  replaced  by  n(n+l).  In  addition,  if  /^(/?)is  replaced  by 

fifR)  =  kR  ^viR)  eq  (15)  can  be  put  in  the  form 

T  dfv  dv  T  T  (  1 V 

Rf^  +  R—+  eR^-  /!  +  -  v  =  0.  (20) 

dR^  dR  I  V  2;  _ 


Thus,  v(R)  is  given  by  a  cylindrical  Bessel  function  of  half  order.  From  this  the  radial 
dependence  can  be  given  as 

(21) 

This  is  actually  a  much  simpler  equation  to  employ  than  the  cylindrical  Bessel  function. 
Because  of  its  frequent  occurrence,  this  is  redefined  as  a  spherical  Bessel  function. 

j„m  =  ^J^^AkR)-  (22) 

Which  leaves  the  final  form  of  the  elementary  wave  function  as 
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fe  =j„ikR)P:(cose)  .  m0.  (23) 

mn  Sin 


This  equation  can  then  be  placed  back  into  the  solution  for  the  E-field  to  give  the  final 


equation  for  the  solution  to  the  spherical  vector  wave  equation. 

^  f  )  w  r  AW 

Rfe 


E(/?,0,0)  = 


1 


dR^ 


+  k  Ry/ 


1 


Rfe 


V  0  J 


kR  dRdO 


)  V 


6  + 


1 


Rfe 


kRsind  dRd(p 


J  V 


<l> 


(24) 
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Appendix  B:  Power  Across  Full  Head  with  Bone  Bolus 


This  appendix  contains  selected  power  profiles  across  the  full  head  with  bone  tissue 
simulation.  The  slices  were  selected  to  show  a  good  representation  of  the  power  profile 
near  the  center  along  with  a  basic  view  of  the  power  near  the  edges.  Note  most  spiking 
occurs  in  the  neck  region  with  a  little  along  the  outside  skin  surface. 
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Z  Location 


Coronal  Cut  Through  the  Head  Model  at  X=38. 


Material  ID  # 


10  20  30  40  50  60  70 

Y  Location 


83 


z  Value  z  Value 


Power  at  X=1 2. 
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Y  Value 
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Z  Value  Z  Value 


Power  at  X=20. 
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Z  Value  Z  Value 


Power  at  X=28. 


Power  at  X=32. 


10  20  30  40  50  60  70 


Y  Value 
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Z  Value  Z  Value 


Power  at  X=34. 
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Z  Value  Z  Value 


Power  at  X=37. 
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Z  Value  z  Value 


Power  at  X=39. 
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Power  at  X=40. 
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Z  Value  Z  Value 


Power  at  X=42. 
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Power  at  X=44. 
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Z  Value  Z  Value 


Power  at  X=48. 


Power  at  X=52. 
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Y  Value 
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z  Value  z  Value 


Power  at  X=56. 
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Z  Value  z  Value 


Power  at  X=64. 


Y  Value 
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ANALYTIC.F  code 


PROGRAM  Analytic 

C  This  program  computes  the  exact  field  distribution  for 
a  sphere 

C  of  tissue  with  a  given  epsr  and  sigma 
INTEGER  COUNT , I , J , K , RADIUS , I 
COMPLEX  JN,EFLDR,EFLDT 
REAL  THETA, PHI , PI , ICNTR, JCNTR, KCNTR 
EXTERNAL  JN,  PN 

C  ***  INITIALIZE  ***  C 
PI=ACOS (-1. ) 

COUNT =0 

CELLSIZE=, 002554 
ICNTR=38 
KCNTR=37 . 5 
JCNTR=38 
RADIUS=37 

OPEN (15 , f ile= ' analytic . dat ’ , status^ ' unknown' ) 

DO  163  1=1,75 
WRITE (6, *)  '1=' ,I 

DO  162  J=l,75 
DO  161  K=l,75 
C 

R=SQRT( (I-ICNTR) **2 . + ( J-JCNTR) **2 . + (K-KCNTR) **2 . ) 

R=R*CELLSIZE 

IF  (R.LT. 0.003)  R=.003 
THETA= ATAN2 ( ( SQRT ((1*1.- ICNTR ) *  *  2 
&  + (J*l. -JCNTR) **2) ) , (K*l. -KCNTR) ) 


PHI=-PI/2 

ELSE 

PHI=ATAN2 ( (J*1,-JCNTR) , (I*1.-ICNTR) ) 
ENDIF 

IF  (PHI.LT.O)  PHI=PHI+2*PI 
IF  (THETA. EQ.O)  THETA= . 00000001 
IF  (PHI. EQ.O)  PHI=. 0000001 


CALL  COMPUTEAP ( R , THETA , EFLDR , EFLDT ) 

C  CALL  CONVERTEF 

(THETA, PHI , EFLDR, EFLDT , EFLDX, EFLDY, EFLDZ ) 

WRITE  (15,*) 

&  sqrt ( (CABS (EFLDR) **2+CABS (EFLDT) **2 ) ) / . 002554 

161  CONTINUE 

162  CONTINUE 

163  CONTINUE 


STOP 

END 


COMPUTEAP  code: 

C  ***  FIELD  AMP  AND  PHASE  COMPUTATION  SUBROUTINE  ******C 

SUBROUTINE  COMPUTEAP  (R, THETA,  EFIELDR, EFIELDT ) 

INTEGER  I,N,IM 

REAL  R,  THETA,  PI , MU, EPS , FREQ,  PN,  OMEGA 
COMPLEX  EFIELDR,  EFIELDT,  AN(5),KR,K,  JN 
EXTERNAL  JN,PN 

PARAMETER ( EPSR= 5 1 , SIGMA=1 . 2 8 ) 
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PI=acos ( -1 . ) 

FREQ=914000000 
0MEGA=2*PI*FREQ 
MU=4*PI*lE-7 
EPS=8 . 854e-12* (EPSR) 

c  K= (0 . , -1. ) *OMEGA*SQRT(MU*EPS) *SQRT ( . 5* (SQRT (1. 

C  &  + (1 .28/ (OMEGA*EPS) ) **2 . ) +1. ) ) 

C  WRITE ( 6 , * )  K 

K= (0 . , -1 . ) * (SQRT ( (0 . , 1. ) *OMEGA*MU* (SIGMA+ (0 . , 1 . ) * 

&  OMEGA*EPS))) 

C  K=IMAG (SQRT (- (OMEGA) **2 . *MU*EPS+ (0,1.) *OMEGA*MU*l . 28 ) ) 
C  WRITE (6,*)  K 

KR=K*R 

C  WRITE (6,*)  KR 

AN(1)= (1.20600400996, .891938522519) 

AN(2)=(0,0) 

AN (3) =(-.19024554426, -.0548692344507) 

AN(4)=(0, 0) 

AN (5)=( .050541055694, -.0191137469021) 

An(2)=(0,0) 

An(4)=(0,0) 

EFIELDR=0 
EFIELDT=0 
DO  500  N=l,3 
I=2*N-1 
IM=I-1 

EFIELDR=EFIELDR+AN{I) * ( (I* (I+l) ) / (KR) ) * 

&  JN(I,KR) *PN(I,COS (THETA) ) 

EFIELDT=EFIELDT+AN (I) * (JN( (IM) ,KR) - 

Sc 

(I/KR) *JN(I,KR) ) * (I/SIN (THETA) ) * (COS (THETA) *PN(I,COS (THETA) ) 
&  PN( (IM) , COS (THETA) ) ) 

500  CONTINUE 
RETURN 
END 
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CONVERTEF  code: 

C  ***  FIELD  CONVERSION  SUBROUTINE  (COMPUTES  FROM  SPHERICAL  TO 
RECT)  ***  C 

SUBROUTINE  CONVERTEF 

(THETA, PHI , EFLDR, EFLDT , EFLDX, EFLDY, EFLDZ ) 

COMPLEX  EFLDR,  EFLDT,  EFLDX,  EFLDY,  EFLDZ 
REAL  THETA, PHI 


EFLDX=SIN (THETA) *COS ( PHI ) *EFLDR+COS (THETA) *COS ( PHI ) *EFLDT 

EFLDY=S IN (THETA) *SIN ( PHI ) *EFLDR+COS (THETA) *SIN(PHI) *EFLDT 
EFLDZ=COS (THETA) * EFLDR- SIN (THETA) *EFLDT 

RETURN 

END' 

C 

'k'k'ie-k'kieir-k'ititiKitir'k'kirir'k'k’k'kitir'kit'kicicieie’kic-kiricie'k'kie'k'kir'k-kitic'k'k'k'kicif-k-kiricifirif^ 


DIRECT.F  Code 


PROGRAM  direct 


C  This  code  is  modified  often  to  produce  plots  of  different 
areas  of  the  exact  source. 

INTEGER  I , J , K , KCNTR , ICNTR , JCNTR , RADIUS 
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REAL  THETA, PI,R, PN, PHI , CELLSIZE, EFXR, EFYR, EFZR 
REAL  EFXI , EFYI , EFZI , MAGEX , MAGEY , MAGEZ 
COMPLEX  EFLDT , EFLDR , JN , EFLDX , EFLDY , EFLDZ 
EXTERNAL  JN,PN 

c  READ  (5,*)  N,  z 

c  WRITE (6,*)  JN(N,Z) 

KCNTR=40 

JCNTR=40 

ICNTR=40 

RADIUS=30 

CELLSIZE=.0032 

R=.0945 

THETA=:ACOS  (-1.  )  12. 

PI=  ACOS (-1. ) 

PHI=PI/2 
DO  100  1=0,  90 
theta= (1/90 . ) *pi 
IF  (R.LT. 0.004)  R=.004 
C  WRITE (6,*)  R, THETA, PHI 

CALL  COMPUTEAP  ( R , THETA , EFLDR , EFLDT ) 

WRITE (6, *) 

( (sqrt (cabs (EFLDT) **2 . +cabs (EFLDR) **2 . ) **2 . ) /2 . ) *1 .28 
C  CALL  CONVERTER ( THETA , PHI , EFLDR , EFLDT , EFLDX , EFLDY , EFLDZ ) 
C  WRITE (6,*)  R,  CABS (EFLDT) 

100  CONTINUE 
STOP 
END 


IDROT.F  Code 


program  idrot 

C  This  program  wioll  rotate  the  ID  matrix  to  make  it 
compatible 
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C  with  the  power  matrix. 

REAL  idl(75,75,75) ,id2(75,75,75) , id3 (75 , 75 , 75 ) 
INTEGER  I,J,K 


OPEN (22, FILE=' S6lD2.dat' , STATUS= ’ OLD ' ) 
OPEN(23,FILE='S6lD3.daf , STATUS =' UNKNOWN ' ) 

DO  101  K=l,75 
DO  102  J=l,75 
DO  103  1=1,75 

READ (22,*)  idl(I,J,K) ,ID2{I,J,K) ,id3 (I,J,K) 

103  CONTINUE 
102  CONTINUE 
101  CONTINUE 

DO  104  1=1,75 
DO  105  J=l,75 
DO  106  K=l,75 

WRITE (23,*)  idl(I, J,K) ,id2 (I, J,K) ,id3 (I, J,K) 
106  CONTINUE 
105  CONTINUE 

104  CONTINUE 
CLOSE(22) 

CLOSE (23) 

STOP 

END 


JN  code: 

C******  Spherical  Bessel  Function  Subroutine  ***********c 
COMPLEX  FUNCTION  JN  (N,Z) 
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INTEGER  N, II 
COMPLEX  Z , J JN (20), JNO 


JN0=CSIN(Z) /Z 
IF  (N.EQ.O)  THEN 
JN=JN0 
RETURN 
END  IF 

JJN(l) =CSIN(Z) /Z**2 .O-CCOS  (Z) /Z 
IF  (N.EQ.l)  THEN 
JN= JJN ( 1 ) 

RETURN 
END  IF 

JJN(2)= (3 . /z**3 .0-1. /Z) *CSIN(Z) - (3 . /Z**2. ) *CCOS (Z) 
IF  (N.EQ.2)  THEN 
JN= JJN ( 2 ) 

RETURN 

ENDIF 

IF  (N.GT.2)  THEN 
DO  100  11=3, N 

JJN(II) = ( (2* (II-l) +1) /Z) *JJN{II-1) -JJN(II-2) 

100  CONTINUE 
JN=JJN(N) 

RETURN 

ENDIF 

RETURN 

END 


MODIFYGEOM.f  Code 

PROGRAM  MODIFYGEOM 


100 


C  THIS  PROGRAM  MODIFIES  THE  GEOMETRY  FILE  TO  SHOW  WHERE 
SOURCES  WERE 
C  SET. 


REAL  R 

REAL  TMP1,TMP2,TMP3,TMP4,TMP5 
INTEGER 

I, J,K, IDONE(92, 92,95) , IDTWO (92 , 92 , 95 ) , IDTHRE ( 92 , 92 , 95 ) 
INTEGER  NX , NY , NZ , ITMPl , ITMP2 
PARAMETER ( ICNTR=48 , JCNTR=48 , KCNTR=4 ) 

PARAMETER (RADIUS=3 7 ) 


OPEN  (22,FILE='SPHERE3.id' , STATUS^ ' OLD ' ) 

OPEN  (23,FILE='SPHERE3.m.id' , STATUS =' UNKNOWN ' ) 
REWIND (23 ) 

READ  (22,*)  NX,NY,NZ 
WRITE(23,*)  NX,NY,NZ 
READ (22,*)  TMPl , TMP2 , TMP3 
WRITE (2 3,*)  TMP1,TMP2,TMP3 
CELLSIZE=TMP1 
READ (22,*)  ITMPl 
WRITE (2 3,*)  ITMPl 
DO  100  1=1,  14 

READ (22,*)  TMPl,  TMP2 ,  TMP3 ,  TMP4 ,  TMP5 
WRITE (2 3,*)  TMPl,  TMP2,  TMP3 ,  TMP4,  TMP5 
100  CONTINUE 

READ (2 2,*)  ITMPl, ITMP2 
WRITE (2 3,*)  ITMPl, ITMP2 

WRITE ( 6 , * )  ICNTR , KCNTR , JCNTR 
DO  101  K=l,  NZ 
WRITE (6,*)  K 
DO  102  J=l,  NY 
DO  103  1=1,  NX 
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READ  (22,  *)  IDONEd,  J,K)  ,  IDTWO  ( I ,  J,  K)  ,  IDTHRE  ( I ,  J ,  K ) 

C  FIRST  CHECK  THE  RADIUS  AT  THE  CENTER  OF  THE  CUBE  OF 
INTEREST 

R=SQRT{( (I*l.+.5)-ICNTR*l.)**2.+( (J*l.+.5)- 
JCNTR*1. ) **2.+ 

&  ( (K*l.+.5) -KCNTR*!. ) **2 . ) 

C  IF  IT  IS  WITHIN  1/2  CELL  OF  TEH  RADIUS  THEN  SET  THE  12  E 
COMPONENTS 

IF( (R.GT, RADIUS-. 5) .AND. { R . LE . RADIUS+ . 5 ) )  THEN 
C  START  BY  SETTING  THE  EX(I,J,K)  COMPONENT 
IDONEd,  J,K)=10 

C  NOW  SET  THE  EXd,J,K+l)  COMPONENT 
IDONEd,  J,K+1)  =10 

C  NOW  SET  THE  EXd,J  +  l,K+l)  COMPONENT 
IDONEd,  J+1,K+1)=0 
C  NOW  SET  THE  EXd,J+l,K)  COMPONENT 
IDONEd,  J+1,K)  =10 
C  NOW  SET  THE  EYd,J,K)  COMPONENT 
IDTWOd,  J,K)=10 

C  NOW  SET  THE  EYd  +  l,J,K)  COMPONENT 
IDTWO (1+1, J,K) =10 

C  NOW  SET  THE  EY ( I+l , J , K+1 )  COMPONENT 
IDTWO (1+1, J,K+1) =0 
C  NOW  SET  THE  EY(I,J,K+1) 

IDTWOd,  J,  K+1)  =10 
C  NOW  SET  THE  EZ(I,J,K)  COMPONENT 
IDTHRE (I,J,K) =10 

C  NOW  SET  THE  EZd  +  l,J,K)  COMPONENT 
IDTHRE (1+1, J,K) =10 

C  NOW  SET  THE  EZ(I+1,J+1,K)  COMPONENT 
IDTHRE (1+1, J+1,K) =10 
C  NOW  SET  THE  EZ(I,J+1,K)  COMPONENT 
IDTHREd,  J+1,K)=10 
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ENDIF 


WRITE(23, *)  IDONE(I,J,K) , IDTWO ( I , J , K) , IDTHRE ( I , J, K) 
103  CONTINUE 
102  CONTINUE 
101  CONTINUE 

CLOSE (22) 

CLOSE (23) 

STOP 

END 


NECKEXTEND.F  Code 

PROGRAM  NECKEXTEND 

C  THIS  PROGRAM  EXTENDS  THE  NECK  DOWN  TO  THE  EDGE  OF  FDTD 

SPACE 


REAL  TMPl,TMP2,TMP3,TMP4,TMP5 
INTEGER 

I, J,K, IDONE(92, 92, 95) , IDTWO (92 , 92 , 95 ) , IDTHRE ( 92 , 92 , 95 ) 
INTEGER  NX,NY,NZ, ITMPl, ITMP2, Z 


OPEN  ( 22 , FILE= ' head5 . id ‘ , STATUS= ' OLD ' ) 

OPEN  (23,FILE='head6.id’ , STATUS UNKNOWN ' ) 
REWIND (23 ) 

READ  (22,*)  NX,NY,NZ 
WRITE(23,*)  NX,NY,NZ 
READ (22,*)  TMPl , TMP2 , TMP3 
WRITE (2 3,*)  TMP1,TMP2,TMP3 
CELLSIZE=TMP1 
READ (22,*)  ITMPl 
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WRITE (2 3,*)  ITMPl 
DO  100  1=1,  14 

READ (22,*)  TMPl,  TMP2 ,  TMP3 ,  TMP4 ,  TMP5 
WRITE (2 3,*)  TMPl,  TMP2 ,  TMP3 ,  TMP4,  TMP5 

100  CONTINUE 

READ (2 2,*)  ITMPl, ITMP2 
WRITE (2 3,*)  ITMPl, ITMP2 

DO  101  K=l,  NZ 
DO  102  J=l,  NY 
DO  103  1=1,  NX 

READ(22,  *)  IDONEd,  J,K)  ,  IDTWO  ( I ,  J ,  K)  ,  IDTHRE  ( I ,  J,  K) 

103  CONTINUE 
102  CONTINUE 

101  CONTINUE 

DO  104  K=l,  17 
WRITE (6,*)  K 
DO  105  J=l,  NY 
DO  106  1=1,  NX 

Z=18-K 

IF  (IDONEd,  J,Z  +  1)  . NE.O. AND. IDONEd,  J,Z  +  1)  .NE.2) 

Sc  IDONEd,  J,Z)=IDONEd,J,Z  +  l) 

IF  (IDTWO (I,J, Z+1) . NE.O. AND. IDTWO(I,J,Z+l) .NE.2) 

Sc  IDTWO  (I,J,Z)=IDTWO(I,J,  Z  +  1) 

IF  (IDTHRE (I,J, Z+1) . NE.O. AND. IDTHRE (I,J, Z+1) .NE.2) 
Sc  IDTHREd,  J,Z)=IDTHRE(I,  J,Z  +  1) 

106  CONTINUE 

105  CONTINUE 

104  CONTINUE 
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DO  107  K=l,  NZ 
DO  108  J=l,  NY 
DO  109  1=1,  NX 

WRiTE(23,  *)  IDONEd,  J,K)  ,  IDTWO  (I ,  J,  K)  ,  IDTHRE  ( I ,  J ,  K) 

109  CONTINUE 
108  CONTINUE 
107  CONTINUE 

STOP 

END 


OPTIMIZE.F  Code 

PROGRAM  Optimize 

C  This  program  computes  the  An  coefficients  for  a  given 

C  eps  and  sigma  at  915  MHZ  for  a  sphere  of  radius  9.45 

cm 


INTEGER  I,N 
COMPLEX  KR, JN,K,AN(5) 

REAL  PI , FREQ , EPS , EPSR, SIGMA, RADIUS , OMEGA, MU , JNP ( 5 ) 
EXTERNAL  JN,  PN 


C  ***  INITIALIZE  ***  C 
PI=ACOS (-1. ) 
RADIUS=.0945 
WRITE (6,*)  'EPSR=?’ 
READ ( 5 , * )  EPSR 
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WRITE (6,*)  'Sigma=?' 

READ (5,*)  sigma 
JNP(1)=.9509 
JNP(3)=.1148 
JNP(5)=.0260 
FREQ=914000000 
0MEGA=2*PI*FREQ 
MU=4*PI*lE-7 
EPS=8 . 854e-12* (EPSR) 

K=  (0 . , -1 . ) * (SQRT ( (0 . , 1. ) *OMEGA*MU* (sigma+ (0 . , 1. ) * 
&  OMEGA*EPS))) 

KR=K*RADIUS 
WRITE(6, *)K,KR 
DO  100  1=1,3 
WRITE (6,*)  'I=’,I 

N=  (1*2) -1 

An(N) =JNP (N) / {JN(N-1,KR) - (N/KR) *JN(N,KR) ) 
WRITE(6,*)  'An(‘,N,‘)=  ■,An(N) 

100  CONTINUE 
STOP 
END 


PN  code: 

C  ***  LEGENDRE  POLYNOMIAL  FUNCTION,  CALCULATES  LEGENDRE 
C  ***  POLYNOMIALS  OF  ORDER  n  (n=0-5) 

REAL  FUNCTION  PN  {N,X) 

INTEGER  N 
REAL  X 
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IF  (N.EQ.O)  PN=1 
IF  (N.EQ.l)  PN=X 
IF  (N.EQ.2)  PN=.5*(3.*X**2.-1. ) 

IF  (N.EQ,3)  PN=0. 5* (5 . *X**3 . -3 . *X) 

IF  (N.EQ.4)  PN= (1. /8 . ) * (35 . *X**4 . -30 . *X**2 . +3 . ) 

IF  (N.EQ.5)  PN= (1. /8 . ) * (63 . *X**5 . -70 , *X**3 . +15 . *X) 

RETURN 

END 


SOURCE.FOR  code  for  the  approximate  source  distribution 

C  FIRST  CHECK  THE  RADIUS  AT  THE  CENTER  OF  THE  CUBE  OF 
INTEREST 

R=SQRT ( ( 1*1 . + . 5-ICNTR) **2 . + ( J*1 . + . 5-JCNTR) **2 . + 
&  (K*l.+.5-KCNTR)**2. ) 

C  IF  IT  IS  WITHIN  1/2  CELL  OF  TEH  RADIUS  THEN  SET  THE  12  E 
COMPONENTS 

IF( (R.GT.RADIUS-.5) .AND. (R.LE.RADIUS+.5) )  THEN 
C  START  BY  SETTING  THE  EX(I,J,K)  COMPONENT 
R=R*CELLSIZE 
IF  ( K . EQ . KCNTR )  THEN 
THETA=ACOS (-1. ) 12. 

ELSE 

THETA=ATAN2 ( (SQRT( ( 1*1 .+. 5-ICNTR) **2. 

Sc  +(J*1.-JCNTR)**2.))  ,  (K*l. -KCNTR)  ) 

ENDIF 

PHI=ATAN2 ( (J*1.-JCNTR) , ( 1*1 .+. 5-ICNTR) ) 

IF  (THETA. EQ.O)  THETA= . 00000001 

IF  (PHI. EQ.O)  PHI=. 0000001 

CALL  COMPUTEAP  (R, THETA, EFLDR, EFLDT ) 
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CALL  CONVERTEF 

(THETA , PHI , EFLDR , EFLDT , EFLDX , EFLDY , EFLDZ ) 

EFXR=REAL ( EFLDX ) 

EFXI=IMAG (EFLDX) 

MAGEX=::CABS  (EFLDX) 

IF  (MAGEX.GT. LIMIT)  THEN 
IF  (EFXR.NE.O.)  THEN 

PHASEX=ATAN2 ( EFXI , EFXR ) 

ELSE 

PHASEX=ACOS (-1. ) /2. 

ENDIF 

NTAU^NINT (T/DT) 

EXS ( I , J , K) =MAGEX/DELX 

&  *FLOAT (MIN(NTAU,NRISE) ) /NRISE*COS (Wl*T+PHASEX) . 

EXS (I, J,K+1) =MAGEX/DELX 

&  *FLOAT (MIN(NTAU,NRISE) ) /NRISE*COS (Wl*T+PHASEX) 

EXS (I, J+1,K+1)=MAGEX/DELX 

&  *FLOAT (MIN(NTAU,NRISE) ) /NRISE*COS (Wl*T+PHASEX) 

EXS ( I , J+1 , K) =MAGEX/DELX 

Sc  *FLOAT  (MIN(NTAU,NRISE)  )  /NRISE*COS  (Wl*T+PHASEX) 

ENDIF 

EFyR=REAL (EFLDY) 

EFYI^IMAG (EFLDY) 

MAGEY=CABS (EFLDY) 

IF  (MAGEY.GT. LIMIT)  THEN 
IF  (EFYR.NE.O.)  THEN 

PHASEY=ATAN2 ( EFYI , EFYR ) 

ELSE 

PHASEY=ACOS(-l.)/2. 

ENDIF 

EYS ( I , J , K ) =MAGEY/DELY 

Sc  *  FLOAT  (MIN  (NTAU,NRISE)  )  /NRISE*COS  (Wl*T+PHASEY) 

EYS  (I  +  l,  J,K)  c=MAGEY/DELY 

Sc  *  FLOAT  (MIN  (NT  AU,NRISE)  )  /NRISE*COS  (Wl*T  +  PHASEY) 

EYS (I+l, J,K+1)=MAGEY/DELY 

Sc  *FLOAT  (MIN(NTAU,NRISE)  ) /NRISE*COS  (Wl*T+PHASEY) 
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phasey=acos (-1 . ) /2 . 

EYS (I, J,K+1) =MAGEY/DELY 

Sc  *FLOAT  (MIN(NTAU,NRISE)  )  /NRISE*COS  (Wl*T+PHASEY) 

ENDIF 

EFZR=REAL (EFLDZ) 

EFZI=IMAG(EFLDZ) 

MAGEZ=CABS (EFLDZ) 

IF  (MAGEZ.GT. LIMIT)  THEN 
IF  (EFZR.NE.O.)  THEN 

PHASEZ=ATAN2 (EFZI, EFZR) 

ELSE 

PHASEZ=ACOS (-1. ) /2 , 

ENDIF 

EZS (I, J,K) =MAGEZ/DELZ 

Sc  *FLOAT  (MIN(NTAU,NRISE)  )  /NRISE*COS  (Wl*T  +  PHASEZ ) 

EZS (I+l, J,K) =MAGEZ/DELZ 

Sc  *  FLOAT  (MIN  (NT  AU,NRISE)  )  /NRISE*COS  (Wl*T+PHASEZ ) 

EZS (I+l, J+1,K)=MAGEZ/DELZ 

&  *FLOAT (MIN(NTAU,NRISE) ) /NRISE*COS (Wl*T+PHASEZ ) 

EZS (I, J+1,K) =MAGEZ/DELZ 

Sc  *FLOAT  (MIN(NTAU,NRISE)  )  /NRISE*COS  (Wl*T  +  PHASEZ ) 

ENDIF 

ENDIF 


SOURCE. FOR  Code  for  an  exact  field  distribution: 


C  FIRST  CHECK  THE  RADIUS  AT  THE  CENTER  OF  THE  CUBE  OF 
C  INTEREST 

R=SQRT ( (I*l.+. 5-ICNTR) **2 .+ (J*l . +. 5-JCNTR) **2 . + 
Sc  (K*l.  +  .5-KCNTR)  **2.  ) 

C  IF  IT  IS  WITHIN  1/2  CELL  OF  TEH  RADIUS  THEN  SET  THE  12  E 
C  COMPONENTS 

IF( (R.GT. RADIUS-. 5) .AND. (R . LE . RADIUS+ . 5 ) )  THEN 
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C  START  BY  SETTING  THE  EX(I,J,K)  COMPONENT 

R=SQRT ( (1*1. +.5-ICNTR) **2 . + (J*l. -JCNTR) **2 . + 

&  (K*l. -KCNTR) **2 . ) *CELLSIZE 

IF  (K.EQ. KCNTR)  THEN 
THETA=ACOS(-l. )/2. 

ELSE 

THETA=ATAN2( (SQRT( ( 1*1 . + . 5 -ICNTR) * *2 . 

Sc  +  (J*l.  -JCNTR)  **2  .  )  )  ,  (K*l .  -KCNTR)  ) 

ENDIF 

PHI=ATAN2 ( (J*l. -JCNTR) , ( 1*1 . + . 5-ICNTR) ) 

IF  (THETA. EQ.O)  THETA= . 00000001 
IF  (PHI. EQ.O)  PHI=. 0000001 
CALL  COMPUTEAP  (R, THETA, EFLDR, EFLDT ) 

CALL  CONVERTER 

( THETA , PHI , EFLDR , EFLDT , EFLDX , EFLDY , EFLDZ ) 
EFXR=REAL (EFLDX) 

EFXI=IMAG (EFLDX) 

MAGEX=CABS (EFLDX) 

IF  (MAGEX.GT. LIMIT)  THEN 
IF  (EFXR.NE.O.)  THEN 

PHASEX=ATAN2 ( EFXI , EFXR ) 

ELSE 

PHASEX=ACOS ( -1 . ) /2 . 

ENDIF 

NTAU=NINT (T/DT) 

EXS ( I , J , K ) =MAGEX/DELX 

&  * FLOAT (MIN ( NT AU,NRISE) ) /NRISE*COS (Wl*T+PHASEX) 

ENDIF 

C  NOW  SET  THE  EX(I,J,K+1)  COMPONENT 

R=SQRT( (1*1, +. 5-ICNTR) **2 .+ (J*l . -JCNTR) **2 .+ 

Sc  (K*l. +1-KCNTR)  **2  .  )  *CELLSIZE 

IF  ( (K+1) .EQ. KCNTR)  THEN 
THETA=ACOS (-1. ) 12. 

ELSE 

THETA=ATAN2 ( (SQRT( ( 1*1 .+. 5-ICNTR) **2 . 

Sc  +(J*1. -JCNTR)  **2.)  )  ,  (K*1.+1-KCNTR)  ) 


no 


ENDIF 

PHI=ATAN2 ( (J*1.-JCNTR) , (1*1 . + . 5-ICNTR) ) 

IF  (THETA. EQ.O)  THETA= . 00000001 
IF  (PHI. EQ.O)  PHI=. 0000001 
■  CALL  COMPUTEAP  (R, THETA, EFLDR, EFLDT ) 

CALL  CONVERTER 

( THETA , PHI , EFLDR , EFLDT , EFLDX , EFLDY , EFLDZ ) 

EFXR=REAL(EFLDX) 

EFXI=IMAG (EFLDX) 

MAGEX=CABS (EFLDX) 

IF  (MAGEX.GT. LIMIT)  THEN 
IF  ( EFXR . NE . 0 . )  THEN 

PHASEX=ATAN2 ( EFXI , EFXR ) 

ELSE 

PHASEX=ACOS (-1. ) /2. 

ENDIF 

NTAU=NINT (T/DT) 

EXS (I, J,K+1) =MAGEX/DELX 

&  *FLOAT (MIN(NTAU,NRISE) ) /NRISE*COS (Wl*T+PHASEX) 

ENDIF 

C  NOW  SET  THE  EX ( I , J+1 , K+1 )  COMPONENT 

R=SQRT( (1*1. +. 5-ICNTR) **2 . + ( J*1 . +1-JCNTR) **2 . + 

&  (K*l. +1-KCNTR) **2 . ) *CELLSIZE 

IF  ( (K+1) . EQ.KCNTR)  THEN 
THETA=ACOS (-1. ) 12. 

ELSE 

THETA=ATAN2 ( (SQRT( ( 1*1 . + . 5 -ICNTR) **2. 

&  + (J*l. +1-JCNTR) **2 . ) ) , (K*l. +1-KCNTR) ) 

ENDIF 

PHI=ATAN2 ( (J*l. +1-JCNTR) , ( 1*1 .+. 5-ICNTR) ) 

IF  (THETA. EQ.O)  THETA= . 00000001 

IF  (PHI. EQ.O)  PHI=. 0000001 

CALL  COMPUTEAP  (R, THETA, EFLDR, EFLDT ) 

CALL  CONVERTER 

( THETA , PHI , EFLDR , EFLDT , EFLDX , EFLDY , EFLDZ ) 

EFXR=REAL (EFLDX) 
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EFXI = IMAG ( EFLDX ) 

MAGEX^CABS (EFLDX) 

IF  (MAGEX.GT. LIMIT)  THEN 
IF  (EFXR.NE.O.)  THEN 

PHASEX=ATAN2 ( EFXI , EFXR ) 

ELSE 

PHASEX=ACOS(-l.)/2. 

ENDIF 

NTAU=NINT (T/DT) 

EXS (I, J+1,K+1) =MAGEX/DELX 

Sc  *  FLOAT  (MIN  (NT  AU,NRISE)  )  /NRISE*COS  (Wl*T  +  PHASEX) 

ENDIF 

C  NOW  SET  THE  EX(I,J+1,K)  COMPONENT 

R=SQRT ( ( 1*1 . + . 5-ICNTR) **2 . + ( J*1 . +1-JCNTR) **2 . + 

Sc  (K*l. -KCNTR)  **2  .  )  *CELLSIZE 

IF  (K.EQ. KCNTR)  THEN 
THETA=ACOS (-1. ) 12. 

ELSE 

THETA=ATAN2 ( (SQRT( ( 1*1 .+. 5-ICNTR) **2. 

Sc  +  (J*1.+1-JCNTR)  **2  .  )  )  ,  (K*l.  -KCNTR)  ) 

ENDIF 

PHI=ATAN2 ( (J*l. +1-JCNTR) , (1*1 .+. 5-ICNTR) ) 

IF  (THETA. EQ.O)  THETA= . 00000001 

IF  (PHI. EQ.O)  PHI=. 0000001 

CALL  COMPUTEAP  ( R , THETA , EFLDR , EFLDT ) 

CALL  CONVERTER 

(THETA, PHI , EFLDR, EFLDT, EFLDX, EFLDY, EFLDZ ) 

EFXR=REAL (EFLDX) 

EFXI = IMAG (EFLDX) 

MAGEX=CABS (EFLDX) 

IF  (MAGEX.GT. LIMIT)  THEN 
IF  (EFXR.NE.O.)  THEN 

PHASEX=ATAN2 ( EFXI , EFXR ) 

ELSE 

PHASEX=ACOS ( -1 . ) / 2 . 

ENDIF 


112 


NTAU=NINT (T/DT) 

EXS ( I , J+1 , K) =MAGEX/DELX 

&  *FLOAT(MIN(NTAU,NRISE) ) /NRISE*COS (Wl*T+PHASEX) 

ENDIF 

C  NOW  SET  THE  EY(I,J,K)  COMPONENT 

R=SQRT( (I*1.-ICNTR) **2 . + ( J*1 . + . 5-JCNTR) **2 . + 

&  (K*l. -KCNTR) **2. ) *CELLSIZE 

IF  (K.EQ. KCNTR)  THEN 
THETA=ACOS (-1. ) /2 . 

ELSE 

THETA=ATAN2 ( (SQRT ( (1*1. -ICNTR) **2 . 

Sc  +  (J*l.  +  .  5-JCNTR)  **2  .  )  )  ,  (K*l .  -KCNTR)  ) 

ENDIF 

PHI=ATAN2 ( (J*l. +. 5-JCNTR) , (1*1. -ICNTR) ) 

IF  (THETA. EQ.O)  THETA= . 00000001 

IF  (PHI. EQ.O)  PHI=. 0000001 

CALL  COMPUTEAP  (R, THETA, EFLDR, EFLDT ) 

CALL  CONVERTER 

( THETA , PHI , EFLDR , EFLDT , EFLDX , EFLDY , EFLDZ ) 

EFYR=REAL(EFLDY) 

EFYI=IMAG (EFLDY) 

MAGEY=CABS (EFLDY) 

IF  (MAGEY.GT. LIMIT)  THEN 
IF  (EFYR.NE.O.)  THEN 

PHASEY=ATAN2 (EFYI,EFYR) , 

ELSE 

PHASEY=ACOS (-1. ) /2 . 

ENDIF 

NTAU=NINT (T/DT) 

EYS ( I , J , K) =MAGEY/DELY 

Sc  *FLOAT(MIN(NTAU,NRISE)  )  /NRISE*COS  (Wl*T  +  PHASEY) 

ENDIF 

C  NOW  SET  THE  EY(I+1,J,K)  COMPONENT 

R=SQRT ( (1*1. +1-ICNTR) **2 . + ( J*1 .+. 5-JCNTR) **2 . + 

Sc  (K*l. -KCNTR)  **2.  )  *CELLSIZE 

IF  (K.EQ. KCNTR)  THEN 
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THETA^ACOS ( - 1 . ) / 2 . 

ELSE 

THETA=ATAN2 ( {SQRT( ( 1*1 , +1-ICNTR) **2. 

&  +  ( J*1 .  +  . 5-JCNTR) **2 . )  ) ,  (K*l . -KCNTR) ) 

ENDIF 

PHI=ATAN2 ( {J*l. +. 5-JCNTR) , (1*1 . +1-ICNTR) ) 

IF  (THETA. EQ.O)  THETA= . 00000001 

IF  (PHI. EQ.O)  PHI=. 0000001 

CALL  COMPUTEAP  (R, THETA, EFLDR, EFLDT) 

CALL  CONVERTER 

( THETA , PHI , EFLDR , EFLDT , EFLDX , EFLDY , EFLDZ ) 

EFYR=REAL (EFLDY) 

EFYI=IMAG (EFLDY) 

MAGEY=CABS (EFLDY) 

IF  (MAGEY.GT. LIMIT)  THEN 
IF  (EFYR.NE.O.)  THEN 

PHASEY=ATAN2 (EFYI, EFYR) 

ELSE 

PHASEY=AC0S(-1.)/2. 

ENDIF 

NTAU=NINT (T/DT) 

EYS (I+l, J,K) =MAGEY/DELY 

&  * FLOAT (MIN (NTAU,NRISE) ) /NRISE*COS (Wl*T+PHASEY) 

ENDIF 

C  NOW  SET  THE  EY ( I+l , J, K+1 )  COMPONENT 

R=SQRT( (I*1.+1-ICNTR) * *2 . + ( J*1 . + . 5 -JCNTR) * *2 . + 

&  (K*l. +1-KCNTR) **2 . ) *CELLSIZE 

IF  ( (K+1) .EQ. KCNTR)  THEN 

THETA=ACOS (-1. ) 12. 

ELSE 

THETA=ATAN2 ( (SQRT( (1*1 . +1-ICNTR) **2. 

&  + (J*l.+. 5-JCNTR) **2. ) ) , (K*l . +1-KCNTR) ) 

ENDIF 

PHI=ATAN2 ( (J*l.+. 5-JCNTR) , ( 1*1 . +1-ICNTR) ) 

IF  (THETA. EQ.O)  THETA= . 00000001 
IF  (PHI. EQ.O)  PHI=. 0000001 
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CALL  COMPUTEAP  (R, THETA, EFLDR , EFLDT ) 

CALL  CONVERTEF 

(THETA, PHI , EFLDR, EFLDT , EFLDX, EFLDY, EFLDZ ) 

EFYR=REAL(EFLDY) 

EFYI=IMAG (EFLDY) 

MAGEY=CABS (EFLDY) 

IF  (MAGEY.GT. LIMIT)  THEN 
IF  (EFYR.NE.O.)  THEN 

PHASEY=ATAN2 ( EFYI , EFYR ) 

ELSE 

PHASEY=ACOS (-1 . ) /2 . 

ENDIF 

NTAU=NINT (T/DT) 

EYS (I+l, J,K+1) =MAGEY/DELY 

&  * FLOAT (MIN (NT AU,NRISE) ) /NRISE*COS (Wl*T+PHASEY) 

ENDIF 

C  NOW  SET  THE  EY(I,J,K+1) 

R=SQRT ( (1*1 . -ICNTR) **2 . + ( J*1 . + . 5-JCNTR) **2 . + 

&  (K*l. +1-KCNTR) **2 . ) *CELLSIZE 

IF  ( (K+1) .EQ.KCNTR)  THEN 
THETA=ACOS (-1. ) 12. 

ELSE 

THETA=ATAN2 ( ( SQRT ((1*1.- ICNTR ) *  *  2 . 

&  + (J*l.+. 5-JCNTR) **2 . ) ) , (K*l .+1-KCNTR) ) 

ENDIF 

PHI=ATAN2 ( (J*l,+. 5-JCNTR) , (1*1. -ICNTR) ) 

IF  (THETA. EQ.O)  THETA= . 00000001 

IF  (PHI. EQ.O)  PHI=. 0000001 

CALL  COMPUTEAP  ( R , THETA , EFLDR , EFLDT ) 

CALL  CONVERTEF 

(THETA, PHI , EFLDR, EFLDT , EFLDX, EFLDY, EFLDZ ) 

EFYR=REAL (EFLDY) 

EFYI=IMAG (EFLDY) 

MAGEY=CABS (EFLDY) 

IF  (MAGEY.GT. LIMIT)  THEN 
IF  (EFYR.NE.O.)  THEN 
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PHASEY=ATAN2 ( EFYI , EFYR ) 

ELSE 

PHASEY=ACOS (-1. ) /2 . 

ENDIF 

NTAU=NINT (T/DT) 

EYS (I, J,K+1) =MAGEY/DELY 

&  *FLOAT {MIN(NTAU,NRISE) ) /NRISE*COS (Wl*T+PHASEY) 

ENDIF 

C  NOW  SET  THE  EZ(I,J,K)  COMPONENT 

R=SQRT ( (1*1. -ICNTR) **2 . + (J*l . -JCNTR) **2 . + 

&  (K*l.+.5-KCNTR)**2. ) *CELLSIZE 

THETA= ATAN2 ( { SQRT ({1*1.- ICNTR ) *  *  2 . 

Sc  +  (J*l. -JCNTR)  **2  .  )  )  ,  (K*l.  +  .5-KCNTR)  ) 

IF  (I.EQ.ICNTR.AND.J.EQ. JCNTR)  THEN 
PHI=-ACOS (-1. ) /2 . 

ELSE 

PHI=ATAN2 ( {J*l. -JCNTR) , (1*1. -ICNTR) ) 

ENDIF 

IF  (THETA. EQ.O)  THETA= . 00000001 

IF  (PHI. EQ.O)  PHI=. 0000001 

CALL  COMPUTEAP  (R, THETA, EFLDR, EFLDT ) 

CALL  CONVERTEF 

( THETA , PHI , EFLDR , EFLDT , EFLDX , EFLDY , EFLDZ ) 

EFZR=REAL (EFLDZ) 

EFZI=IMAG (EFLDZ) 

MAGEZ=CABS (EFLDZ) 

IF  (MAGEZ.GT. LIMIT)  THEN 
IF  (EFZR.NE.O.)  THEN 

■PHASEZ=ATAN2 (EFZI, EFZR) 

ELSE 

PHASEZ=ACOS (-1. ) /2. 

ENDIF 

NTAU=NINT(T/DT) 

EZS ( I , J , K) =MAGEZ/DELZ 

&  *FLOAT (MIN(NTAU,NRISE) ) /NRISE*COS (Wl*T+PHASEZ ) 

ENDIF 
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C  NOW  SET  THE  EZ(I+1,J,K)  COMPONENT 

R=SQRT ( (1*1. +1-ICNTR) **2 . + ( J*1 . -JCNTR) **2 . + 

&  (K*l.+.5-KCNTR)**2.)*CELLSIZE 

THETA=ATAN2( (SQRT( (1*1 . +1-ICNTR) **2. 

Sc  +  (J*l. -JCNTR)  **2  .  )  )  ,  (K*l.  +  .5-KCNTR)  ) 

IF  (  (I  +  l)  .EQ.ICNTR. AND. J.EQ. JCNTR)  THEN 
PHI=-ACOS (-1. ) /2. 

ELSE 

PHI=ATAN2 ( (J*l. -JCNTR) , (1*1 . +1. -ICNTR) ) 

ENDIF 

IF  (THETA. EQ.O)  THETA= . 00000001 

IF  (PHI. EQ.O)  PHI=. 0000001 

CALL  COMPUTEAP  (R, THETA, EFLDR, EFLDT ) 

CALL  CONVERTEF 

(THETA, PHI , EFLDR, EFLDT, EFLDX, EFLDY, EFLDZ ) 

EFZR=REAL (EFLDZ) 

EFZI=IMAG (EFLDZ) 

MAGEZ=CABS (EFLDZ) 

IF  (MAGEZ.GT. LIMIT)  THEN 
IF  (EFZR.NE.O.)  THEN 

PHASEZ=ATAN2 (EFZI,EFZR) 

ELSE 

•  PHASEZ=ACOS (-1. ) /2. 

ENDIF 

NTAU=NINT (T/DT) 

EZS (I+l, J,K) =MAGEZ/DELZ 

Sc  *FLOAT(MIN(NTAU,NRISE)  )  /NRISE*COS  (Wl*T+PHASEZ ) 

ENDIF 

C  NOW  SET  THE.  EZ (I+l, J+1,K)  COMPONENT 

R=SQRT( (I*1.+1-ICNTR)**2.+(J*1.+1.-JCNTR)**2.+ 
Sc  (K*l.  +  .5-KCNTR)**2.)*CELLSIZE 

THETA=ATAN2 ( (SQRT( (1*1 . +1 . -ICNTR) **2. 

Sc  +  (J*l.+1. -JCNTR)  **2.  )  )  ,  (K*l.  +  .5-KCNTR)  ) 

IF  (  (I  +  l)  .EQ. ICNTR. AND.  (J+1)  .EQ. JCNTR)  THEN 
PHI=-ACOS (-1. ) /2 . 

ELSE 
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PHI=ATAN2 ( (J*1.+1-JCNTR) , (1*1 . +1 . -ICNTR) ) 

ENDIF 

IF  (THETA. EQ.O)  THETA= . 00000001 
-  IF  (PHI. EQ.O)  PHI=. 0000001 

CALL  COMPUTEAP  (R, THETA, EFLDR, EFLDT) 

CALL  CONVERTEF 

( THETA , PHI , EFLDR , EFLDT , EFLDX , EFLDY , EFLDZ ) 

EFZR=REAL ( EFLDZ ) 

EFZI=IMAG (EFLDZ) 

MAGEZ^CABS (EFLDZ) 

IF  (MAGEZ.GT. LIMIT)  THEN 
IF  (EFZR.NE.O.)  THEN 

PHASEZ=ATAN2 (EFZI, EFZR) 

ELSE 

PHASEZ=ACOS(-l.)/2. 

ENDIF 

NTAU=NINT (T/DT) 

EZS (I+l, J+1,K)=MAGEZ/DELZ 

&  *FLOAT (MIN(NTAU,NRISE) ) /NRISE*COS (Wl*T+PHASEZ ) 

ENDIF 

C  NOW  SET  THE  EZ(I,J+1,K)  COMPONENT 

R=SQRT( (1*1. -ICNTR) **2.+(J*l.+l.-JCNTR)**2.+ 

Sc  (K*l.  +  .5-KCNTR)**2.)*CELLSIZE 

THETA=ATAN2 ( (  SQRT ((1*1.- ICNTR ) *  *  2 . 

Sc  +  (J*l.  +1.  -JCNTR)  **2  .  )  )  ,  (K*l .  +  .  5-KCNTR)  ) 

IF  (I. EQ. ICNTR. AND. (J+1) .EQ. JCNTR)  THEN 
PHI=-ACOS (-1. ) /2 . 

ELSE 

PHI=ATAN2 ( (J*1.+1-JCNTR) , (1*1. -ICNTR) ) 

ENDIF 

IF  (THETA. EQ.O)  THETA= . 00000001 

IF  (PHI. EQ.O)  PHI=. 0000001 

CALL  COMPUTEAP  ( R , THETA , EFLDR , EFLDT ) 

CALL  CONVERTEF 

(THETA, PHI , EFLDR, EFLDT , EFLDX, EFLDY, EFLDZ ) 

EFZR=REAL (EFLDZ) 
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EFZI=IMAG (EFLDZ) 

MAGEZ=CABS (EFLDZ) 

IF  (MAGEZ.GT. LIMIT)  THEN 
IF  (EFZR.NE.O.)  THEN 

PHASEZ=ATAN2 (EFZI,EFZR) 

ELSE 

PHASEZ=ACOS (-1 . ) /2 . 

ENDIF 

NTAU=NINT (T/DT) 

EZS (I, J+1,K) =MAGEZ/DELZ 

&  * FLOAT (MIN (NT AU,NRISE) ) /NRISE*COS (Wl*T+PHASEZ ) 

ENDIF 


ENDIF 


WRITEEF  code: 

C  ***  SUBROUTINE  TO  WRITE  OUT  FIELDS  TO  A  FILE  ******** *c 

SUBROUTINE  WRITEEF  ( EFLDX, EFLDY , EFLDZ , I , J , K , COUNT ) 

COMPLEX  EFLDX, EFLDY, EFLDZ 
REAL 

EFXR, EFXI , EFYR, EFYI , EFZR, EFZI , TESTX, TESTY, TESTZ , LIMIT 
REAL  PHASEX, PHASEY, PHASEZ 
INTEGER  I,J,K,COUNT 

LIMIT=.01 

OPEN (UNIT  =  2  5 , FILE= " SPHERESOURCE . DAT " , STATUS= " OLD " , ACCESS= ' APP 
END'  ) 

EFXR=REAL (EFLDX) 
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EFYR=:REAL  (EFLDY) 

EFZR=REAL (EFLDZ) 

EFXI=IMAG (EFLDX) 

EEYI=IMAG (EFLDY) 

EFZI=IMAG (EFLDZ) 

TESTX=CABS (EFLDX) *1000 
TESTY^^CABS  (EFLDY)  *1000 
TESTZ=CABS (EFLDZ) *1000 
IF  (TESTX.GT. LIMIT)  THEN 
COUNT=COUNT+l 
IF  (EFXR.NE.O)  THEN 

PHASEX=ATAN2 (EFXI,EFXR) 

ELSE 

PHASEX= (ACOS (-1. ) /2) 

END  IF 

WRITE (25, 80)  I , J, K, 0 , TESTX, PHASEX 
END  IF 

IF  ( TESTY. GT. LIMIT)  THEN 
COUNT =COUNT+l 

IF  (EFYR.NE.O)  THEN 
PHASEY=ATAN2 (EFYI, EFYR) 

ELSE  - 

PHASEY=AC0S ( - 1 . ) / 2 

END  IF 

WRITE(25,80)  I, J,K,1,TESTY,PHASEY 
ENDIF 

IF  (TESTZ.GT. LIMIT)  THEN 
COUNT=COUNT+l 

IF  (EFZR.NE.O)  THEN 
PHASEZ=ATAN2 (EFZI,EFZR) 

ELSE 

PHASEZ=ACOS (-1. ) /2 

ENDIF 

WRITE (25, 80)  I, J,K,2,TESTZ,PHASEZ 
ENDIF 
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80  F0RMAT(I3,I3,I3,I2,E14.6,E14.6) 

CLOSE(25) 

RETURN 

END 

Q  *********************★******★****************************(;> 
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Appendix  D:  Matlab  Programs 
AVG.M  Code 

function  y^avg (file, file2) 

%  This  function  computes  the  average  material  characteritics 
%  throughout  the  head/bolus  structure.  File  is  the  ID  file. 

%  File2  is  a  file  containing  the  electrical  characteristics 
%  of  each  material. 

%  Initialize 
top=size ( f ile) 
count (1) =0; 
count (2 ) =0 ; 
count (3)=0; 
y(l,l)=0; 
y (l,2)=0; 
y(l,3)=0; 
y(2,l)=0; 
y(2,2)=0; 
y(2,3)=0; 

%  Run  through  the  locations  and  the  x,y  and  z  components 
%  keeping  a  running  tally  of  permitivity  and  permeability, 
for  i=l:top(l,l) 
for  j=l:top(l,2) 

%  If  the  location  isn't  free  space  (Is  part  of  the  head/bolus 
%  structure) . 

if  file(i,j)>l, 

y(l,j)=y(l,j) +f ile2 (file (i, j ) -1, 1) ; 
y(2,j)=y(2,j) +f ile2 (file (i, j ) -1, 2)  ; 
count ( j ) =count ( j  )  +1 ; 
end 
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end 

end 


%  Compute  the  average  across  all  three  components 
Y(l,4)=(y(l,l)+Y(l,2)+y(l,3))/( count ( 1 ) +count ( 2 ) +count (3)); 
Y(2,4)=(y(2,l)+Y(2,2)+Y(2,3) )/ (count ( 1 ) +count (2 ) +count (3 ) ) ; 

%  Compute  the  x,y,  and  z  averages 
for  j=l:top(l,2) 

Y(l,j)=Y(l,j) /count  (j); 

Y(2,j)=Y(2,j)/ count ( j ) ; 
end 


AXY.M  Code 

function  Y=axY ( f ile , idf ile, level, s tat ) 

%  this  function  will  analyze  an  xy  level  of  an  output  file 
%  First  plot  the  ids  for  the  level.  file  is  the  power  file, 
%  ID  file  is  a  file  of  the  id's,  level  is  the  level  to  be 
%  analyzed,  stat  is  a  string  file  that  describes  the  files 
%  (for  display) 
idpxy (idf ile, level, stat) ; 

%now  plot  the  power  picture 
pwrpxy (file, level, stat) ; 

%  Now  plot  the  normalized  power, 
meshxyn ( file, level, [ -60  45], stat); 
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AYZ.M  Code 


function  y=ayz (file, idf ile, level, stat) 

%  this  function  will  analyze  an  yz  level  of  an  output  file 
%  First  plot  the  ids  for  the  level,  file  is  the  power  file, 
%  ID  file  is  a  file  of  the  id's,  level  is  the  level  to  be 
%  analyzed,  stat  is  a  string  file  that  describes  the  files 
%  (for  display) 
idpyz (idfile, level, stat) ; 

%now  plot  the  power  picture 
pwrpyz (file, level, stat) ; 

%  Now  plot  the  normalized  power, 
meshyzn ( file, level, [ -60  45] , stat) ; 


CLEARSOURCE.M  Code 

function  y=clearsource (meshfile,x) 

%  This  function  clears  everything  from  the  source  out  for  a 
%  particular  cut  of  the  sphere,  meshfile  is  a  file  of  a 
%  single  layer  of  the  power  file.  X  is  the  level  of  the 
%  layer. 

y=f ile; 

for  i=l:75 
for  j=l:75 

%  Compute  the  radius  at  the  point 

r=sqrt( ( i-37 . 5 ) "2+ ( j -37 . 5 ) ^2+ (x-37 . 5 ) ^2 ) ; 

%  if  the  point  is  outside  th  esource,  replace  it  with  an  NaN 
if  r>36 . 5000001, 
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Y (i,  j ) =NaN; 
end 
end 
end 


DIFF.M  Code 

function  y=DIFF (file, level, type) 

%  this  function  computes  the  id  changes  across  two  level 
%  changes.  file  is  th  eidfile  of  interes,  level  is  the  level 
%  of  interes,  and  type  is  a  string  to  describe  the  scenario 
%  (for  display).  Insert  (:,1,2  or  3)  after  the  TEST*  files 
%  to  plot  the  cahnges  of  only  the  x,y  or  z  componenets. 

TESTl=meshxy (file, level-1) ; 

TEST2=meshxy (file, level) ; 

TEST3=meshxy ( f ile, level+1)  ; 

figure 

pcolor (TEST1-TEST2 ) 
colormap (gray) 
colorbar 

ylabel ( ' Y  Value ' ) 
xlabel ( ' X  Value ' ) 

title ([' Changes  in  ',type,'  From  Z=' ,num2str (level-1) , ‘  to 
Z= ' , num2str ( level ),’.']) 


figure 

pcolor (TEST3 -TESTl ) 
colormap (gray) 
colorbar 

ylabel ( ' Y  Value ' ) 
xlabel ( ' X  Value ' ) 
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title ([' changes  in  ■,type,'  From  Z= ' , num2str ( level+1 ) , '  to 
Z= ' ,num2str (level) 


ERRYZ  Code: 

function  y=erryz (filel, file2, level) 

%  This  function  computes  the  error  along  one  YZ  slice  of  two 
%  files.  Filel  is  the  “correct"  fiel  and  file2  is  the  file 
%which  is  being  compared,  level  is  the  level  of  interest. 
y= zeros (75,75) ; 

tmpl=clearsource (meshyz (filel, level) , level) ; 
tmp2=clearsource (meshyz (file2, level) , level) ; 

for  i=l;75 
for  j=l:75 

if  tmpl (i , j ) ~=NaN, 

y  (i,  j  )  =  (  (tmpl  (i,  j  )  -tmp2  (i,  j  )  )  /  tmpl  (i,  j  )  )  *100  ; 
else 

y ( i , j ) =NaN ; 
end 
end 
end 

mesh  (y) 


ERRXY.M  Code 

function  y=errxy (filel, f ile2 , level) 
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%  This  function  computes  the  error  along  one  XY  slice  of  two 
%  files.  Filel  is  the  “correct"  fiel  and  file2  is  the  file 
%which  is  being  compared,  level  is  the  level  of  interest. 
y= zeros (75,75) ; 

tmpl=clearsource (meshxy ( filel, level) , level) ; 
tmp2=clearsource (meshxy ( file2 , level) , level) ; 

for  i=l:75 
for  j=l:75 

if  tmpl ( i , j ) ~=NaN, 

y  (i,  j  )  =  (  (tmpl  (i,  j  )  -tmp2  (i,  j  )  )  /  tmpl  (i,  j  )  )  *100  ; 
else 

y (i, j ) =NaN; 
end 
end 
end 

mesh  (y) 


EXPAND.M  Code 

function  y = expand (FILE, THRSH) 

%  This  function  takes  a  file  and  sets  any  value  above  a 
%  threshold  to  30  in  order  to  make  the  hot  points  standout. 
y=FILE; 
for  i=l;75 
for  j=l:75 
if  y (i, j ) >THRSH, 
y(i, j)=30; 
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end 

end 

end 


GETAXIS.M  Code 


function  y=getaxis ( f ile) ; 

%  This  function  computes  the  maximum  and  minimum  of  a  power 
%  file  in  order  to  get  values  for  the  axes  of  a  movie  plot. 
SIZE=size ( file) ; 

ZMAX=1; 

ZMIN=0; 

for  i=l;SIZE(l,l) 

f  ilel=sqrt  ( f  ile  ( i ,  1)  ''2  +  f  ile  (i,  2 )  ''2  +  f  ile  (i ,  3  )  ^2 )  ; 
if  filel>ZMAX, 

ZMAX=f ilel; 
end 

if  filel<ZMIN, 

ZMIN=f ilel; 
end 
end 

y=[ZMIN  ZMAX] ; 


IDPXY.M  Code 

function  y=idpxy (idf ile, level, stat) 
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%  This  function  plots  a  picture  of  the  specified  XY  plane  of 
%  an  ID  file,  stat  is  a  string  describing  the  IDfile. 
figure 

idxy=clearsource (meshxy (idfile, level) , level) ; 

idxy=idxy ; 

colormap (hsv) ; 

colorbar 

pcolor ( idxy) ; 

xlabel  'X  Value' 

ylabel  'Y  Value' 

title  (['Material  id''s  across  an  XY  cut  at 
Z= ' , num2str (level) ,'.']) 
text (15 , -7 ,[' Results  for  ',stat]); 


IDPYZ.M  Code 

function  y=idpyz ( idfile, level, stat) 

%  This  function  plots  a  picture  of  the  specified  level  of  an 
%  idfile.  stat  is  a  string  description  of  the  idfile. 
figure 

idyz=clearsource (meshyz (idfile, level) , level) ; 

idyz=idyz*75/6 ; 

colormap (hsv) ; 

colorbar; 

pcolor (idyz ) ; 

xlabel  ‘ Z  Value ' 

ylabel  'Y  Value' 

title  (['Material  id''s  across  an  YZ  cut  at 
X= ' , num2str (level) ,'.']) 
text (15 , -7 ,[' Results  for  ',stat]); 
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MESHXY.M  Code 


function  y=meshxy ( file, pin) 

%  This  file  creates  a  mesh  in  the  specified  xy  plane 
%  It  also  tests  the  file  and  plots  the  total  field  if 
%  multiple  fields  are  given  in  the  input. 

counts (pin) ; 
test=size ( file) ; 
for  i=l:75 
for  j=l:75 
if  test(2)==l, 
y (i, j ) =f ile (count) ; 
end 

if  test(2)==2, 

y  (i,  j  )  =sqrt  (file  (count ,  1)  ''2  +  f ile  (count ,2)^2)  / 2-, 
end 

if  test(2)==3, 

y  (i,  j  )  =sqrt  (file  (count ,  1)  ''2  +  f ile  (count,  2  )  ■^2  + 
file (count ,3)^2)  ; 
end 

count=count+75 ; 
end 
end 


MESHXYN.M  Code 

function  y=meshxyn ( f ilel , x, view, TXT ) 

%  This  function  creates  a  mesh  of  the  specified  yz  plane 
%  and  plots  the  mesh  fully  formatted. 
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y=nrml (meshxy ( f ilel , x) , vlu (file 1,38, 38, 38) ) ; 
figure 

y^clearsource (y,x) ; 
colormap (hsv) ; 
mesh (y, view) 
xlabel  ' X  Value ' 
ylabel  'Y  Value' 

zlabel  ' |E|^2  normalized  to  1  at  center' 
title  {['E  Field  Distribution  Across  XY  ' , num2str (x) , ' 
Plane . ' ] ) 

msg= [' Results  for  ',TXT]; 
text  ( -40 , 35 , 0 , msg) 


MESHYZ.M  CODE 

function  y=meshyz (file, pin) 

%  This  file  creates  a  mesh  in  the  specified  yz  plane 
%  It  also  tests  the  file  and  plots  the  total  field  if 
%  multiple  fields  are  given  in  the  input. 

count = (pln-1) *75*75+1; 
test^size ( file) ; 
for  j=l:75 
for  k=l:75 
if  test(2)==l, 
y ( j , k) =f ile (count) ; 
end 

if  test(2)==2, 

y(j,k)=file (count , 1) +f ile (count , 2 ) ; 
end 

if  test(2)==3. 


131 


y  ( j  ,  k)  =sqrt  (file  (count ,  1)  ^2  + file  (count ,  2  )  ^^2  + 
f  ile  ( count ,  3  )  ^^2  )  /2  ; 

end 

count=count+l; 

end 

end 


MESHYZN.M  Code 

function  y=meshyzn ( file, x, view, TXT) 

%  This  function  creates  a  mesh  of  the  specified  yz  plane 
%  and  plots  the  mesh  fully  formatted. 

y= (meshy z (file,x) ) ; 

y=clearsource (y , x) ; 

figure; 

colormap (hsv) ; 
mesh (y , view) ; 
shading  faceted; 
xlabel  'Z  Value’ 
ylabel  'Y  Value' 
zlabel  ' |E|^2' 

title  (['E  Field  Distribution  Across  YZ  Plane  at 
X= ' , num2str (x) ] ) 
msg= [' Results  for  ',TXT]; 
text  ( -69 , 35 , -17 ,msg) 
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MMOVIE.M  Code 


function  y=mmovie ( f ile, name, A) 

%  This  function  creates  a  movie  through  the  xy  slices  of  a 
%  file.  This  allows  cycling  through  the  levels  quickly  to 
%  spot  hotspots  easily. 

Y=moviein (75 ) ; 

for  i=l:75 

j=i 

pcolor (meshxy (file, j ) )  ; 
caxis (A); 
colorbar; 
colormap ( jet) ; 

title  ([name,'  at  Z= ', num2str ( j ),'.'] ) 
xlabel  'Y  Location' 
ylabel  'X  Location' 
y ( : , i) =get frame; 
end 

movie  (y ) 


PLOTCY.M  Code 


function  y=plotcy ( f ilel, f ile2 ,x, Z) 

%  This  fnction  plots  the  power  along  a  y  line  along  with  the 
%  material  id's  along  the  line.Filel  is  the  idfile,  file2  is 
%  the  power  file,  and  X  and  Z  define  the  line  of  interest. 

tmpl=meshyzid ( f ilel ( : , 1) ,x) ; 
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tmp2=meshyzid ( f ilel ( : , 2 ) , x) ; 
tmpS^^meshyzid  (  f  ilel  (  :  ,  3  )  ,x)  ; 
tmp4=meshyz (f ile2 ,x) ; 
k=x2  ; 

1=1:75;' 

plot (I, tmpl (:,k),‘x',I, tmp2 (:,k), 'o', I, tmp3 ( : , k) , ' + ' , 
I , tmp4 ( ; , k) , ' - ' ) 
xlabel  'Y  Value' 
ylabel  'Power  (W/m^2)' 

tmp3= [' Total  Power  with  material  lD''s  at  X= ' ] ; 
title  ( [tmp3 ,num2str (x) , '  Z= ' ,num2str (x2 ) , ' . ‘ ) 


PLOTCZ.M  Code 

function  y=plotcz (filel, file2,x,x2) 

%  This  function  plots  the  power  along  a  Z  line  along  with  the 
%  material  id's  along  that  line  different  files,  filel  is  the 
%  id  file,  file2  is  the  power  file  of  interest.  X  and  Y 
%  define  the  line  of  interest. 


tmpl=meshyzid ( filel ( : , 1) ,x) ; 
tmp2=meshyzid ( filel ( : , 2 ) ,x) ; 
tmp3=meshyzid ( filel (:,3),x); 
tmp4=meshyz (file2,x) ; 
k=x2 ; 

1=1:75; 

plot (I, tmpl (k, : ) , ' o ' , I, tmp2 (k, tmp2 (k, : ) , 
'x ' , I, tmp4 (k, ) 
xlabel  ' Z  Value ' 
ylabel  ' Power ' 

tmp3=['Plot  of  Power  and  Material  ID's  at  X= ' ] ; 
title  ( [tmp3 ,num2str (x) , '  Y= ', num2str (x2 ),'.'] ) 
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PLOTX.M  Code 


function  Y=plotx ( file,x,x2 , type) ; 

%  This  function  plots  a  power  file  along  a  y  line.  X  is  the 
%  Y  value,  X2  is  the  Z  location  of  the  line.  type  is  a 
%  string  giving  the  style  of  line  to  draw. 
tmp=meshxy (file,x2) ; 
k=x; 

plot  (tmp( : ,k) ,type) ; 
xlabel  ' Y  Value ' 
ylabel  'E-field  strength' 

title  (['E  field  along  line  X= ' , num2str (x) , ' , 

Z= ' ,num2str (x2) 


PLOTY.M  Code 

function  y=ploty (file, x,x2, type) ; 

%  This  function  plots  the  power  file  along  a  y  line.  X  is 
%  the  S  value,  X2  is  the  Z  value  and  type  is  the  type  of  line 
%  to  plot  the  line  as. 
tmp=meshyz (file,x) ; 
k=x2  ; 

plot  (tmp ( ; ,k) , type) ; 
xlabel  ' Y  Value ' 
ylabel  'E-field  strength' 

title  (['E  field  along  line  X= ' , num2str (x) , ' , 

Z= ' , num2str (x2 ),’.']) 
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PLOTZ.M  Code 


function  Y=plotz (file, x,x2, type) ; 

%  This  function  plots  the  power  file  along  a  y  line.  X  is 
%  the  X  value,  X2  is  the  Y  value  and  type  is  the  type  of  line 
%  to  plot  the  line  as. 
tmp=meshyz (file,x) ; 
k=x2  ; 

plot  (tmp (k, : ) , type) ; 
xlabel  ‘Y  Value' 
ylabel  'E-field  strength' 

title  (['E  field  along  line  X= ' , num2str (x) , ' , 

Z= ' , num2str (x2 ),'.’]) 


POWERAVG.M  Code 

function  y=poweravg (pf ile, idf ile, matpar ) 

%  This  function  computes  the  average  power  at  the  center  of 
%  the  head  as  well  as  across  the  entire  head.  It  also 
%  computes  the  power  over  three  different  thresholds,  pfile 
%  is  a  file  of  the  E-field  over  the  head,  idfile  is  the  id's 
%  of  the  head  mesh,  matpar  is  a  file  of  th  ematerial 
%  characteristics  for  each  id. 

Tmp ( : , 1) =matpar (round (idf ile ( : , 1) ) +1) ; 

Tmp ( : , 2 ) =matpar ( idf ile { : , 2 ) +1) ; 

Tmp ( : , 3 ) =matpar (idfile ( :  ,  3 ) +1)  ; 

power =pf ile. *Tmp; 

powerl=sqrt  (power  ( :  ,  1)  .  ^2+power  ( : ,  2 )  .  ''2+power  ( :  ,  3  )  .  ''2 )  ; 
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power l=powerl/vlu ( power 1, 38,38,38) ; 


%  Compute  the  average  power  at  the  center. 

Tmp=0 ; 
count =0 ; 
for  i=37:39 
for  j=37:39 
for  k=37 : 39 

Tmp=Tmp+vlu  (powerl ,  i ,  j  ,  k)  ; 
count=count +1 ; 
end 
end 
end 

centeravg=Tmp/ count 
Y (1, 1) =centeravg; 

%  Now  compute  the  average  power  across  the  entire  head 
hpower=powerl (~isnan (powerl) ) ; 
to talavg=mean ( hpower ) 
y ( 1 , 2 ) =totalavg; 

y(l,3)=y(l)/y(2)  ; 

%  Now  compute  the  power  spikes 
SIZE=size (hpower) 
numspikes=0; 
num9spikes=0 ; 
num8spikes=0 ; 
count =0 ; 

for  i=l:SIZE(l, 1) ; 
count=count+l ; 
if  hpower (count )>=. 8 , 
num8spikes=num8spikes+l; 

if  hpower ( count )>=. 9 , 
num9spikes=num9spikes+l; 
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if  hpower (count ) >=1, 
numspikes=numspikes+l  ; 
end 
end 
end 
end 

numspikes 

y  ( 1 , 4 ) ^numspikes  ; 

num9 spikes 

y (1, 5) =num9 spikes ; 

numSspikes 

y ( 1 , 6 ) =num8  spikes ; 


PWRPXY.M  Code 


function  y=pwrpxy (file, level, st at ) 

%  This  functions  plots  a  picture  of  the  power  across  the 
%  level  specified,  stat  is  a  string  decripter  to  be  printed 
%  on  the  plot. 

pxy=clearsource (nrml (meshxy (file, level) , 
vlu (file, 38,38,38) ) , level) ; 

figure 

colormap (hsv) ; 
pooler (pxy ) ; 
colorbar ; 
brighten ( . 75 )  ; 
xlabel  'X  Value' 
ylabel  'Y  Value' 

title  (['E  field  magnitude  across  an  XY  cut  at 
Z= ' , num2str (level) ,'.']) 
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text (15 , -7 ,[' Results  for  ■,stat]); 
text (74.5,77, ' E  field  Mag ' ) ; 


PWRPYZ.M  Code 

function  Y=pwrpY2 (file, level, st at ) 

%  This  function  plots  apicture  of  the  power  across  the  YZ 
%  plane  specified,  stat  is  a  string  descriptor  to  be  used  in 
%  the  title  of  the  plot. 

PYz=clearsource (meshYz (file, level) , level) ; 
figure 

colormap (hsv) ; 
pcolor  (pYz;)  ; 
brighten ( . 75 )  ; 
colorbar ; 
xlabel  'Z  Value' 

Ylabel  'Y  Value' 

title  (['E  field  magnitude  across  an  YZ  cut  at 
X= ' , num2str (level ),'.']) 
text (15, -7, [ 'Results  for  ',stat]); 
text (74 . 5 , 77  ,  ' E  field  Mag'); 


SHRINK.M  Code 

function  Y=shrink ( f ile) 

%  This  function  shrinks  an  idfile  from  the  92X92X95  grid  down 
%  to  a  75X75X75  grid  to  match  the  E-field  files. 


count =0 ; 
counts=0 ; 
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Y= zeros (421875,3)  ; 
for  k=l:92 
disp (k) 
for  j=l:95 
for  i=l:92 
counts=counts+l; 

if  (i>10  &  i<86  &  j>10  &  j<86  &  k>10  &  k<86), 
count=count+l ; 
y (count, ; ) =file (counts, : ) ; 
end 
end 
end 
end 
end 


TERR.M  Code 

function  y=terr (filel, file2) 

%  this  function  calculates  teh  error  between  two  files  across 
%  all  of  space. 
y= zeros (211873,1) ; 
count =0 ; 
countl=0 ; 
for  i=l:75 
i 

for  j=l:75 
for  k=l:75 
count=count+l; 

if  sqrt ( (i-38)^2+{j-38)^2+(k-38)^2)<37, 
count l=countl+l; 
y (countl) =abs ( ( filel (count) - 

f ile2 (count) ) /filel (count ) ) *100 ; 
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end 

end 

end 

end 


VLU.M  Code 


function  y=vlu ( f ile , x, y, z ) 

%  This  function  finds  the  value  of  the  file  at  the  specified 
%  location. 
test=size ( file) ; 

count = (x-l)*5625+(y-l)*75+(z-l)+l; 
if  test(2)==l, 
y=f ile (count , 1) ; 
end 

if  test(2)==3, 

y=sqrt  (file  (count,  1)  ^2  +  f  ile  (count ,  2)^2  +  file  (count ,  3  )  ''2 )  ; 
end 
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